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THESIS DISCLAIMER 






The reader is cautioned that computer programs developed 
in this research may not have been exercised for all cases of 
interest. While every effort has been made, within the time 
available, to ensure that the programs are free of 
computational and logic errors, they cannot be considered 
completely validated. Any application of these programs 
without additional verification is at the risk of the user. 
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TABLE I DEFINITION OF TERMS 



KE 

PE 

6 

F 

U 

r. 

R 

d 

W 

U 

W 

"e 

w 

r 

D 

« 

r 

c 

9L 

^## 

'<1 

«* 



Kinetic Energy 
Potential Energy 

Large motion generalized coordinate 
Large motion generalized force, applied moment 
column vector for small motion generalized 
coordinates . 

local position vector (L subscript refers to load) 

global position vector (L subscript refers to load) 

deformation vector of transverse deflections 

mass density of the link (L subscript refers to load) 

transformation matrix (a function of 0) 

derivative of the transformation matrix wrt 0 

residual accelerations term 

deformation transformation matrix 

shape function matrix derived in Appendix C 

special derivative of the shape function matrix 

stiffness matrix 

gravitational acceleration vector 

volume integral with a kernel for dependent matrices 
where ## represents an index number from Table A 
length integral (along the boom) 
volume integral where the # represents an index 
number # from Table A. 2 
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TABLE II 



DEFINITIONS OF INTEGRAL EXPRESSIONS 



^11 

^12 

^21 

^22 



(w'^W) = Ju (W*^ 

(W'^W) = Ju (w"^ 

(W*^W) = Ju (W*^ 

(W*^W) = Ju (W*^ 

= J c r dx 



W) £. dv 



W) dv 



W) r dv 



W) $ dv 



(the kernel would be an 
indicated function of W) 



= Ju dv 
"4 ' f“L ^ 



Hg (d'^w'^WDj^) 



/ 



UlcJ<d’'m^MDk) c, 



dv. 



(given the tip load as a 
Point load, and the 
kernel as function of 
D and W) 
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I. INTRODUCTION 



A. BACKGROUND 

The U.S. Navy has declared it's intention to enhance its 
global, real-time monitoring of the ocean characteristics and 
to improve its ability to forecast oceanographic parameters 
by means of a dedicated space system [Ref. 1]. The system 
under consideration would improve the battle group 
commander's ability to screen his forces from detection by 
enemy submarines/ be more aware of severe weather conditions 
in the operating theater and improve forecasters meteorology 
predictions . 

The satellite will use a variety of special sensors and 
operate from an 820 km polar orbit. A number of physical 
considerations for this type of mission embody the foreground 
of modern structural and control engineering research. 

One such item considered necessary to achieve the Navy's 
desired goals is a Low Frequency Microwave Radiometer (LFMR). 
It represents a unique aspect of spacecraft engineering and 
control. The LFMR would provide sea surface temperatures 
accurate to within a half of a degree Kelvin with a surface 
resolution capability of 25 km. The LFMR would require a 
large (approximately 5 meter diameter) antenna with precise 
pointing knowledge; its field of view would scan at either 31 
or 16 scans per minute. The width of scan would be at 
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approximately 1400 km centered on the ground track, in front 
of or behind the satellite. The incident beam would make a 
53.1 degree angle with the local vertical at the surface, see 
Figure 1. [Ref. 2:p. 3] 

The large ground coverage is provided by the 
eccentrically mounted antenna. A schematic drawing of the 
proposed LFMR sub-system required to support the antenna is 
presented in Figure 2. The boom sizes and off-axis angles 
vary from contractor to contractor, those represented in 
Figure 2 are not considered the only feasible choices. 

B. MOTIVATION 

Currently the Navy .has halted development on the Navy 
Remote Ocean Sensing System (N-ROSS), which was to be the 
vehicle for the LFMR. Cost considerations and the 
technological development (simulations and modeling ) of the 
LFMR boom assembly and control scheme have been suggested as 
primary concerns. 

This investigation is motivated by the limited number of 
techniques currently employed for the dynamic analysis of 
large flexible space structures. We will present the 
development and computer simulation of an analysis for a 
flexible spacecraft boom, rotating in three dimensions. 
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C. LITERATURE REVIEW 



The approach to modeling the kinematics and dynamics of 
the rotating flexible space antenna are not unlike those 
techniques employed in the field of robotics necessary for 
the design and control of flexible manipulators. In the case 
of flexible manipulator models, the focus is to develope the 
relationship between the large, rigid-body motion (large 
motion) and the small motions due to the structure's 
flexibility . 

The most popular class of these reduction methods for 
dynamic analysis of structures appears to employ the modal 
superposition technique [Refs. 3-7]. In this type of 
analysis a determination must be made to limit the number of 
contributing mode shapes and then combine those that will 
contribute in a way which yields a system of equations of 
motion. 

Book [Refs. 4,5] utilizes this modal technique to model a 
system with flexible kinematics. He utilizes recursive 
Lagrangian kinematics and large and small motion coordinates. 
The resulting equations of motion are an accurate model 
describing large and small motion and the coupling between 
them. However the result, a set of coupled, non-linear 
second order, ordinary differential equations, proved 
computationally intensive, and therefore expensive to solve. 

Sunada and Dubowsky [Ref. 8] on the other hand expressed 
the kinematics and dynamics of the manipulator in terms of 
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"small perturbations about some known nominal motion of the 
system". The effect of the coupling between the large and 
small motion is incomplete when the large motion is assumed 
nominal and the small motion superimposed upon it. 

Chang's development of the Equivalent Rigid Link System 
ERLS [Ref. 9], describes the large motion kinematics of the 
system by the ERLS and the small motion deformations relative 
to it. The application of finite element techniques and 
Lagrangian dynamics, produces two sets of coupled, non- 
linear, ordinary differential equations of motion. The 
application of the ERLS formulates these sets of equations; 
one set for large motions and one set for small motions. The 
set of large motion equations are non-linear in both the 
large and small motion variables. The set of small motion 
equations are linear in the small motion variable and non- 
linear in the large motion variables. 

The characteristics of these equation sets permit a 
solution by numerical integration with a variety of 
techniques possible. The integration technique developed by 
Chang, the Sequential Integration Method, is straightforward 
and relatively easy to implement. Chang's model offered a 
complete representation of the dynamics of flexible 
manipulators as well as one that could be efficiently solved 
using the Sequential Integration Method. 

Simo and Vu-Quoc [Ref. 10] present an approach which 
appears well suited for design but the choice of generalized 
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coordinates may require a complicated observer to relate the 
computational prediction of the generalized coordinates with 
observable reference quantities of the controller 
environment . 

Previous work at the Naval Postgraduate School utilizing 
the ERLS technique validated the procedure in a planar 
application [Refs. 11,12]. Both Petroka and Gannon applied 
Chang's ERLS utilizing separate mode shape response matrices. 
The agreement between the assumed modal responses and the 
experimentally observed response justifies this first effort 
at applying the ERLS to model the dynamic response of 
flexible booms in the spacecraft applications. 

D. OUTLINE 

The nonlinear dynamic response of a spacecraft boom 
structure member is predicted by the numerical integration of 
the equations of motion. In this work, the spatial finite 
element discretization of the boom and the application of an 
assumed polynomial mode shape provides an approximate 
solution. 

The modeling process has been applied through the ERLS 
utilizing Lagrangian Dynamics to obtain the system equations 
of motion for numerical integration. 

Chapter Two of this work reviews Chang's development of 
the ERLS and in Chapter Three the ERLS is applied to a 
rotating flexible boom and a variety of simulation results 
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are presented. 


The computer simulation code and methodology 


are introduced 


in Chapter Four. The results of this work 


illustrate both 


the need for controlling the flexibility and 


the versatility 


of the ERLS solution technique. A review of 



this work and suggested areas for development is presented in 
Chapter Five. 
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II. DYNAMIC M ODEL FORMULAILQILEQR.A ..SEACECRAFT 
SINGLE-ELEMENT FLEXIBLE BOOM 



A. KINEMATIC RELATIONSHIPS OF THE EQUIVALENT RIGID LINK 

SYSTEM 

The introduction and development of the Equivalent Rigid 
Link System, separates the kinematics of the flexible body 
into large and small motions. The large motions are modeled 
with the characteristics of rigid bodies while the small 
motions are modeled by finite element techniques and defined 
relative to the position of the equivalent rigid bodies. 

The small motions resulting from loading flexure are 
measured in the ERLS relative to the local coordinate system 
for each equivalent rigid body. The absolute position of a 
point on a deformed body is obtained by coordinate 
transformations which utilize the characteristics of joint 
variables as in the methods of Hartenberg and Denavit [Ref. 
131. Lagrangian dynamics are used to obtain the equations of 
motion and therefore time derivatives of the absolute 
position will be necessary in the development of the 
Lagrangian expression. 

The ERLS model accounts for the change in location 
experienced by each point along a flexible member. Discrete 
points or nodes are chosen and their relative position 
determined while locations of intermediate points are 
assumed according to a mathematical function. The choice of 
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nodes and mathematical function determines the theoretical 



shape of the flexed member. 

The Finite Element Method (FEM) is used to descretize the 
flexible member assigning the nodes and enforcing 
compatibility between elements. The choice of shape function 
and node position must be consistent with the accuracy 
desired. The ERLS will accommodate closed form natural mode 
shape functions or any order of polynomial. In this work for 
example a cubic polynomial shape function was chosen to 
represent the bending displacement of the flexible boom. The 
nodal positions chosen here are the end points of the 
flexible member. The selection of a • polynomial shape 
function and end point nodal positions are consistent with 
previous works which have experimentally validated their 
selection. [Refs. 9,11,12]. 

The kinematics of the ERLS separates .the global motions 
into the sums of large and small motions. The large motion 
is described by the equivalent rigid link and the small 
motion due to deformations of the link are described relative 
to the rigid link. 

The characteristics of rigid body dynamics is essentially 
classical dynamics. In this application the Hartenberg- 
Denavit 4x4 transformation matrix is used to completely 
describe the positions and velocities of each point in the 
system. Appendix A develops these kinematic relationships in 
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a brief presentation, a complete review of the kinematics of 
the ERLS is available [Ref 9:pp. 17-501. 

B. KINETIC RELATIONSHIPS OF THE EQUIVALENT RIGID LINK 

SYSTEM 

The formulation of the Lagrange equations to derive the 
equations of motions is straightforward and systematic. This 
familiar approach is helpful in the complex analysis of this 
type. 

The selection of generalized coordinate systems is 
accomplished indirectly by the separation of the global 
motion into large and small motions. Formalizing this choice 
assigns the large motion generalized coordinate as theta (9), 
which represents a rotation and the small motion generalized 
coordinates (U) a column vector composed of the translations 
and rotations at the boom tip. 

The translations from the nodal position on the ERLS to 
the deformed position on the actual flexible member (x,y,z) 
will become the first three small motion generalized 
coordinates. The slopes at the node or the rotational angles 
of flexible member about the axis perpendicular to the planes 
of displacement of the flexible member (U') are included as 
the last three small motion generalized coordinates. 

The transformation from the ERLS to the deformed boom 
position is accomplished by the mode shape matrix. The 
derivation of the mode shape matrix is briefly presented in 
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Appendix B. The choice of nodal position and mode function 
determines the accuracy of the simulation response. 

In this representation a cubic polynomial function was 
chosen as the shape function for bending and a linear 
function for the extension and torsion. This choice provides 
the mathematical representation of our first three small 
motion generalized coordinates, and simplifies the required 
mathematics to obtain the functional representation of the 
remaining three small motion generalized coordinates. 

The system kinetic energy is the sum of the kinetic 
energies of each member of the system, this would include the 
flexible boom and the tip load (if any). The system 
potential energy results from the elastic strain energy of 
the flexible member, and the application of gravitational 
body force (if any). Generalized forces represent the 
application of applied forces and damping forces. Appendix C 
includes the development of the equations of motion from the 
Lagrange Equations for the case of end point nodal positions 
and a polynomial shape function for six small motion 
generalized coordinates. 

The manipulations of the Lagrange equations to obtain the 
equations of motion entail extensive effort, mathematical 
dexterity, and fierce determination but the reward is a set 
of two non-linear, coupled, second-order, ordinary 
differential equations. Though both sets are coupled, one 
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describes the large motions and the other describes the small 



motions, they are presented below. 



M e’ + M U' = F (2-1) 

qq qn q 

M e" + M U'+ G U + K U = £ (2-2) 

nq nn n n n 



where : 
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n 



1x1 effective inertia matrix for large motions 
1x6 coupled effective inertia matrix of the 
small motion on the large motion 
1x1 load vector for the large motion 
6x1 coupled effective inertia matrix of the 
large motion on the small motion 

6x6 effective inertia matrix for small motions 
6x6 gyroscopic matrix 
6x6 stiffness matrix 

6x1 load vector for the small motion 
generalized coordinate of the large motion 
6x1 generalized coordinate vector for small 
motion deformations 



C. NUMERICAL INTEGRATION OF THE EQUATIONS OF MOTIONS 

The application of numerous numerical integration schemes 
is possible in the solution of the equations of motion. The 
primary consideration in the selection of the integration 
technique is the determination of the time step necessary to 
integrate the equations of motion and obtain numerical 
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stability. In this work several integrations techniques were 
tested i.e.: Runga-Kuta 5th order (both fixed and variable 
time step), Adams-Moul ton (variable time step), and Gear. 
While all of these techniques achieved numerical stability 
the computational time required to obtain adequate simulation 
time was excessive. 

The formulation of the equations of motion in a manner 
similar to equations (C-36) and (C-37), permits the 
application of the Sequential Integration Method [Ref. 9:pp 
79-91J. In this technique an implicit integration scheme is 
employed to integrate the small motion equation, and an 
explicit technique is used to integrate the large motion 
equation. 

The implicit technique guarantees numerical stability in 
the case of linear equations and hence an iterative technique 
can be avoided. The results of the small motion integration 
can be provided to the large motion equation which can then 
be integrated explicitly in a predictor-corrector fashion. 
The choice of time step in the case of the explicit 
integration can be governed by solution accuracy vice small 
motion frequency response. The decided advantage to this 
approach is the computational time savings. 
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III. APPLICATION OF THE ERLS TO A SPACECRAFT 
SINGLE ELEMENT FLEXIBLE BOOM 



A. THE LOW FREQUENCY MICROWAVE RADIOMETER (LFMR) 

The low frequency microwave radiometer (LFMR) which would 
be based on a space platform of the N-ROSS type is composed 
of a reflector, an upper reflector boom, a lower reflector 
boom, and a spacecraft boom. These booms would all be hinged 
together and owing to the requirements of space systems, 
would be developed from light weight materials. The dynamic 
modeling of such a structure has been the primary focus of 
the work presented here, other applications of the ERLS in 
space technology do however exist. 

The orientations of the microwave receiving dish on board 
requires* precise pointing accuracy which must be well 
controlled for the transmission/reception of signals. The 
accuracy of the pointing mechanism is subject to the position 
error determinations which in turn are vulnerable to the 
accelerations of orbital adjustment, antenna spinning and the 
inherent light weight materials of construction. For these 
reasons the control of vibrations is essential. 

The vibrational response or deformation resulting from 
the application of these forces must be observable to the 
controller. The choice of the small motion coordinates 
previously outlined would make these deformations available 
during start up, regular operation, and orbital station 
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Since these values 



keeping adjustments. Since these values are directly 
measurable this choice could simplify the controller observer 
requirements while ensuring the mission specifications are 
met . 



B. SIMULATION OBJECTIVES 

In this work the multi-bodied boom has been simplified to 
a single member. The deformations resulting from the 
structural flexibility, have been assumed small and hence 
small angle assumptions have been implemented in the model. 
The application of body forces and internal forces have been 
limited to gravity, and a driving torque. No damping effects 
have been implemented, nor contributions to deformations from 
thermoelastic effects upon entering or leaving the eclipse. 
The incorporation of these concerns is possible in a similar 
model but beyond the scope of the work presented here. 

The exact material of the boom would likely be a graphite 
epoxy composite with specific material properties, however 
the simulation code utilizes the material considerations of 
pure aluminum to permit the intuitive interpretations 
associated with a more familiar structure material. Specific 
material considerations can be easily changed. 

The ultimate goal of this work was to complete the first 
stages of the simulation modeling. To develope the solution 
techniques and the basic simulation code (Appendices A-D) 
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which can later be expanded to eventually include multi- 
elements in the satellite model. 



C. DESCRIPTION OF THE COMPUTER SIMULATION CODE 

The IBM simulation language DSL (Dynamic Simulation 
Language) was originally chosen to model the equations of 
motions. The variety of integration schemes and FORTRAN 
compatibility were the over riding reasons for this choice. 
Gannon's work in two dimensions indicated a preferred 
integration technique for the* Adams-Moulton method [Ref. 
12;p. 25], additionally DSL proved capable of simulating 
Gannon's system model. 

While integration methods were not specifically the focus 
of this work, as indicated previously, the choice of 
integration technique became dependant upon the time step 
requirements for numerical stability. 

The simulation code was developed with modular FORTRAN 
subroutines then attached to a short initialization program. 
From the three dimensional development in Appendix C the 
Lagrange equations can be formulated as equations (2-1) and 
(2-2) and then expanded in the form of equations (C-36) and 
(C-37). Where the formulation in this manner lends itself to 
the Sequential Integration Method. 

The resulting column vector; 



(e 



u u u u u 
X y z XX yy zz 



(3-1) 
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represents the large motion acceleration 6 and the rate of 
local deformation in position and slope. Where the x y and z 
subscripts represent the local coordinates and the double x y 
and z subscripts represent the rotations about these axis. 

Each subroutine of the simulation code was developed to 
specifically determine a portion of the single element three 
dimensional analysis and due to the nature of modular 
programing the expansion of this code to a multi-body problem 
can be accommodated. The logic flow of the coding 
(Figure 15) is provided to assist in the future development 
of multi-body simulations. 

The simulation code can be viewed from three levels; 
LEVEL I (an overview) separates the code into the primary 
focal points of INITIALIZATION, FORMULATION, INTEGRATION, 
REFORMULATION, and, OUTPUT. The reformulation is a technique 
to incorporate the inertia coupling between the large and 
small motion equations. 

LEVEL II can be viewed as the FORTRAN formulation of 
coefficients matrices, and sequential integration to solve 
for ©, I, 6*, U, U, and u" . 

LEVEL III is a listing of the FORTRAN utilities sub- 
programs required for the matrix formulation and manipulation 
in LEVEL II. Each of the units in LEVEL II may access the 
routines in LEVEL III. 
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D. SIMULATION RESULTS 



The computer simulation code was run with a variety of 
initial conditions these conditions were chosen to present 
the aspects of the boom flexibility and the impact this 
flexibility would have on the requirements for end point 
control. The simulation output has been presented in local 
displacement values, where these relative values represent 
the conditions at the boom tip, the point mass location. 

The simulation results presented in Appendix E are 
considered typical of results available from the simulation 
code. They represent the single large motion generalized 
coordinate and each of the six small motion generalized 
coordinates. The intent of these figures is to provide the a 
physical feeling for the plant behavior. 

The results. Figures 3-7, reflect the slowing effects on 
the large motion generalized coordinates from the point mass. 
The small motion effects. Figures 8-14, and the coupling 
between the large and small motion coordinate systems, is 
also clearly observable. Additionally the results indicate 
the effects of boom elevation angle and corresponding 
responses at four elevations, zero through 6QO. 
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IV. SUMMARY 



A. CONCLUSIONS 

This investigation is motivated by the needs of the 
United States Navy to be able to accurately predict the 
Meterological phenomenon within an operations theater. To 
this end the techniques currently employed for the dynamic 
analysis of large flexible space structures must be improved 
and expanded. Ne have presented the development and computer 
simulation of an analysis for a flexible spacecraft boom, 
rotating in three dimensions with flexible modes in bending, 
extension and torsion. 

The nonlinear dynamic responses of the structure member 
are predicted by the numerical integration of the equations 
of motion. In this work, the spatial finite element 
discretization of the boom structure and the application of 
an assumed third order modal response provide an approximate 
solution to the equations of motion, one which lends itself 
to the problems associated with controller requirements. 

The modeling process has been applied through the 
Equivalent Rigid Link System, utilizing Lagrangian Dynamics 
to obtain the system equations of motion for numerical 
integration. The effectiveness of the integration schemes 
available for numerical integration of this type of problem 
have been briefly presented. 
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The justification for the 



ultimate selection of a 



Sequential Integration Method was discussed. Areas of the 
computer code which must be developed in order to extend the 
model to multiple bodies and eventually to the controller 
design and application were also presented. 

B. RECOMMENDATIONS 

Areas which remain to be investigated utilizing the 
simulation code include: 

1. The significance of torsion and extension in the 
total system. 

2. The significance of the way in which torque is 
applied to the boom, for example; impulse, ramp or step 
wise changes. 

3. The significance of the material properties of the 
model . . 

4. The significance of the rotational velocity. 

The single largest concept remaining to be addressed by 
this boom model is the additional members to be joined to the 
driven boom. The addition of these members requires 

additional application of the ERLS modeling scheme. 
Implementation requires compatibility for both members at 
the connecting joints. 

The long-term goal of this research is the design and 
control of a large flexible antenna-boom assembly for a space 
based meterological data acquisition system. Simulation 
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studies need to continue and their scope should be expanded 
to multi-body analysis. 

The control system requirements must be investigated and 
incorporated into the antenna system model. Finally a 
physical model must be developed to validate the simulation 
results . 
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APPENDIX A 



DERIVATION OF THE KINEMATICS 
REQUIRED IN THE EQUATIONS OF MOTION 



In Appendix C we presented the derivation of the 
equations of motion for the single-element, flexible boom. 
In the derivation of these equations we represent the 
position of any point in the deformed member globally by 
applying either a coordinate transformation matrix W or a 
deformation transformation matrix D and in the determination 
of the small motions the application of both transformations 
is required. 

In order to more accurately represent the coordinate 
transformation matrix W, we will present the Hartenberg- 
Denavit transformation matrix A, and demonstrate the 
technique for obtaining the matrix W and it's derivatives. 



A 



1 

a COS 0 
a SIN e 
d 



0 

COS 0 
SIN 0 
0 



0 

-SIN0 COSa 
COS0 COSa 
SIN a 



0 

SIN0 SINa 
-COS0 SINa 
COS a 



(A-1) 

The theory and derivation of the Har tenberg-Denavit 
matrix formulations are defined in several of the references 
cited in Chapter One this work [Refs. 9,11-13]. The intent 
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here is to present the application of this technique. 
Variables variables used in the H-D transformation matrix are 
called joint variables and these are defined below. 

6 = the rotation of the link about the local z axis 
a = the length of the flexible link 
d = the z axis displacement between consecutive 
coordinate systems 

a E the angular rotation between consecutive local z 
axis 



In the application of the H-D matrix, the rotation angle 
0, of the member, is always about the local z axis. To 
obtain this orientation we have defined two sets of local 
coordinates which will require two coordinate 
transformations. The first transformation is a rotational 
one, we transform the global z axis to the local y axis, 
(this is a rotation of « radian about the x axis). The 
transformation matrix A1 represented in equation (A-2), 
provides this rotational transformation. 



A 



1 



1 0 

0 COS 0 

0 SIN 0 

0 0 



0 

-SIN0 COSa 
COS0 COS« 
SINa 



0 

SIN0 SIN« 
-COS0 SINa 
COS0 



(A-2) 
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The second local coordinate system is defined at the end 
of the link, and therefore we have included the effects of 
translation. We have represented this as A 2 it is defined 
by equation (A-3). 



1 

a COS 8 
a SIN 8 
0 



0 0 
COS 8 -SIN8 COSa 

SIN 8 C0S8 COSa 

0 SIN a 



0 

SIN8 SINa 
-COS8 SINa 
COS 8 

(A-3) 



Having defined the two H-D matrices required in the 
simulation, we can represent the coordinate transformation 
matrix W. 



W 




(A-4) 



The definition of W lends itself toward the application 
of the chain rule in the determination of the derivatives of 
W. By the choice of the H-D transformation matrix we can 
easily obtain the derivatives of each A matrix as follows: 




(A-5) 



Where the derivative is taken with respect to the 
generalized coordinate, in this application q would be the 
angle theta and Q is defined by equation (A-6), 
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Q 



(A-6) 



0 0 0 0 

0 0-10 

0 10 0 

0 0 0 0 



If we applied the unique conditions of the spacecraft 
rotating boom we can set the change in Q 2 as zero and 

eliminate the concern for order in multiplication of the 
three matrices which compose the derivative of the W matrix. 



and : 




QA1A2 



W. 




0 



(A-7) 



(A-8) 



For the time derivatives we obtain (with 6 constant); 



A 



dW 

dt 




90 _ 9W 

at ae. 



■ 

e. 



W - ^t W QA^0j^A2 QA^0^0^A2 + ^2 



(A-9) 



(A-10) 



In expression (A-10) we can consolidate the product of 

each group of terms on the extreme right hand side and 

■ 2 

rename them the residual acceleration 0^^ and the large 
motion acceleration W 0*. In following derivations 
subscripts for 0 have been omitted since there is only one 
large motion generalized coordinate capable of motion, the 
derivatives of the other would be zero. 
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Defining the 


deformation transformation 


matrix as 


equation 


(A-11) 


will 


accommodate three 


dimensional 


de formation. 








D = 


1 

“k=i 

'^K=2 

*^K=3 


0 

1 

“k=6 

~“k=5 


0 0 

-"k=6 “k=5 

1 -"k=4 

“k=4 1 


(A-11) 



The partial derivative of D with respect to each of the 
small motion generalized coordinates results in the same 
matrix as the partial of the time derivative of D with 
respect to the time derivative of the small motion 
generalized coordinates, equation (A-12). In the following 
development the subscript K refers to the particular small 
motion generalized coordinate of interest. 






PD 

PUk 



9D 




(A-12) 
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APPENDIX B 



DERIVATION OF THE MODE SHAPE RESPONSE MATRIX ♦ 

The derivation of the mode shape matrix has been 
presented previously [Refs. 9,11,12]. The derivation 
presented here is abbreviated and included only as a brief 
description of the techniques involved, for complete 
derivation the reader is encouraged to review the previous 
works . 

Representing the relative position changes due to 
deformations by d. we have assumed small deformations, no 
warping after torsion and all deformations u, v, w, and 0 are 
only functions of x (along the boom), the yz plane remains 
plane after torsion and deformations only vary along the x 
direction. 





" 0 














u 


y 


+ z 0 

y 








1 


V 


^ ®x 








(B-1) 




w + 


y ®x 










Applying 


simple 


beam 


theory. 


we have 


assumed that the 


rotary effect and 


the 


lateral 


shear 


deflections are 


negl igible . 
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0 



+ y 0 

X 

We have chosen a polynomial form to represent the shape 
function for the displacement of the single element flexible 
boom as shown below. 



"l 




^0 


+ 


aj^x 


(B- 3 ) 


•^2 


= 






2 3 

bj^x + b2X + b^x 


(B- 4 ) 




= 


*=0 


+ 


2 3 

C^X + C2X + C3X 


(B- 5 ) 




= 


^1 






(B-6) 




= 






2 

2b2X + Sb^x 


(B- 7 ) 




= 


^1 




2 

2C2X + 3C3X 


(B-8) 



We can now apply the boundary conditions of zero 
displacement and slope at the base where x is equal to minus 
the boom length, L, and the sign is consistent with the 
convention for the local coordinate system. 

Substituting the boundary conditions into the shape 
functions gives: 
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0 



(B-9) 





-L) 


— 


^0 


+ 


aj^x = 




-L) 






+ 


bj^x + 




-L) 


= 




+ 


Cj^X + 




-L) 


= 


^1 








-L) 


S 


‘^l 


+ 


2 b2X 




-L) 


= 


*=1 


+ 


2 C 2 X 



'=2*^ * '’ 3 *^ 

+ 3 
+ 3 



0 

0 



(B-10) 

(B-11) 

(B-12) 

(B-13) 

(B-14) 



Uj^(O) 

U2(0) 

U3<0) 

u^(0) 

«5(0) 

U6(0) 




(B-15) 

(B-16) 

(B-17) 

(B-18) 

(B-19) 

(B-20) 



The coefficients a 2 » ^ 2 ' *^3' ^2' ^ 3 ' 

solved for by substitution and arranged in the shape matrix 
representation as; 







0 

0 

0 



'45 




(B-21) 



where: R = — ^ x, the variable of integration, L the length 
over which the integration will take place (the length of the 
flexible member). 
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1 + R 



^21 = 

332 = - 2R^ - 3R^ + 1 

a^3 = - 2R^ - 3R^ + 1 

33 ^ = - 2 (1 + R) 

= - Y (1 + R) 

a^g = LR^ + 2LR^ + LR 

a 3 g = LR^ + 2LR^ + LR 

For the potential energy it will be necessary to develop 
the expression for u|, u^, u^, u^,u^, u^, accomplished by 
differentiating the shape function equations (B-9) - (B-20) 
and then solving the simultaneous equations a second time to 
produce the modified shape function matrix which is presented 
below : 



r 



Cj^^ 0 0 0 0 

“ “ =23 “ =25 

0 C32 000 

000 c.. 0 

44 



0 

0 



'36 



0 



(B-22) 
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1 








"T" 




®23 




12R 


6__ 












6R 


4 


^25 




“IT ■ 


L 


®32 


= 


°23 




®36 




^25 




*=44 


= 


^11 
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APPENDIX C 



DERIVATION OF THE EQUATIONS OF MOTION FOR THE 
SINGLE-ELEMENT, FLEXIBLE BOOM 



In derivations of this type additional definitions are 
made to simplify the terms involved in the Lagrange 
equations. A table of terms, for both variable names and 
integral definitions (Tables I, and II respectively) have 
been included for this purpose. 

In the case of a rotating flexible boom the large and 
small motions can be separated by application of the ERLS. 
The large motion generalized coordinate is defined as d and 
the small motion generalized coordinates as U. In this 
application U represents a column vector of six entries. Ue 
can imagine this vector of generalized coordinates as 
composed of two sets. The first set ul, u2, and u3 represents 
the translational displacement of the end point in the local 
X y and z directions when deformed. The second set u4, u5, 
and u6 represents the rotation or slope of the flexed member 
at the deformed end. In accordance with this description the 
vector of generalized coordinates can be defined by equation 
(C-1) . 



T 




*^1 “2 *^3 “4 “5 *^6 



] 



(C-1) 
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These generalized coordinates will be used in the 
application of the Lagrange equations to develope the 
equations of motion for both the large and small coordinate 
systems. The Lagrange equations are defined in equations 
(C-2) and (C-3). 




Small motion generalized forces are zero in this example,, 
this is consistent for a system without damping. Utilizing 
the ERLS we can define the global position of the. end point 
in terms of the local position vector (r) and a deformation 
vector (d) equation (C-4). Similarly the global position of 
the tip load is obtained by transforming a local position 



vector for the 


load, vector 


(rL), 


first 


by 


a deformation 


transformation 


matrix D 


and 


then 


by 


a coordinate 


transformation 


matrix W. 











E. = W (r + d) 



(C-4) 



WDt:^ 



(C-5) 



In order to completely represent the coordinate 
transformation matrix we have presented the Hartenberg- 
Denavit transformation matrix A, and demonstrate the 
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technique for obtaining the matrix W and it's derivatives in 
Appendix A. Thus we continue with the development of the 
equations of motions, not encumbering the derivation of these 
equations with the derivations of others. Also presented in 
Appendix A is the definition of the deformation 
transformation matrix D. 

Continuing; we present the derivative with respect to 
time for the global positions in order to obtain the kinetic 
energy expressions. 

R = W(r + d) + Wd (C-6) 

= WDr^+ WOr^ (C-7) ’ 

The kinetic and potential energy expressions equations 
(C-8) and (C-9) follow. In this development we have 
separated the large and small motion energy contributions and 
will present them separately. 

For the large motion: 

KE = j [ Jw R dv + Ju ^ dvj (C-8) 

link load 

volume volume 

PE = J Ju'^r'^C ru dx - J |i g dv - J rJ g dVj^ (C-9) 

link link load 

length volume volume 
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Now redefining d in terms of the shape matrix $ (again 
this derivation is presented in Appendix B to preserve the 
continuity here). We can rewrite equations (C-9) as: 

KE = Ju[j[(r+*U)'^w'^W(rL+$U) + (*U)'^W^W*U] + ( r+*U ) *^w'^W*U] dv 

+ J)m[i[(Dr^'^w'^W(Dr^ + ( Dr^'^w'^WC Dr^ ] + (Dr^'^w'^W(Dr^ ] dVj_^ 

(C-10) 

Continuing in the development we first express the 

derivative of the kinetic energy with respect to the time 

■ 

rate of change of the large motion coordinate, 0, equation 
(C-11), and then obtain the time rate of change of this 
expression, equation (C-12). 

= fu[(r+*U) Vw(r+*U) + (r+*U)Vw*5] dv + 

90 j •- 0 0 j 

JuL[(Dr^'^W^W(Dr^ + (Dr^’^WgW(Dr^) ] dVj^ (C-11) 



dfC^] ' Ju[(^U)'^wJw(r+^U) + (r+$U)'^(WQW+WQW‘) (r+*U) + 
(r+*U)'^(wJw) (*0) + ($5 )'^(WqW) (*U) + 
(r+*U)'^(WQW+wJw) (§U) + (r+*U)'^(WQW) (♦()')] dv + 

Jy^[(Dr^'^WQW(D£^ + (Dr^'^(WQW + Wg5')(Dr^ + 

(Dr^^wjfl(6£i^ + (6r^'^WQW(6r^ + ( Di:^'^WqW( 6'r^ + 



(Dr^'^(WQW + WQW)(Dri^] dv^^ 



(C-12) 
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Finally we express the derivative of the kinetic energy 
with respect to the large motion position as 0 . 

=Ju[|[(E.+*U)'^(WqW + w'^W^Xr+^U) + (^U)'^(WqW + w'^W^H^U)] 

+ (r+*U)'^(WQW+fl’’wQ) ($U)] dv + Jw[j[(Dr^’’(® 0 W + ®'^We>(Dr^ + 

(Sii^TWqW + w'^W^HDc^] f (Dl^TWqW + w'^WQ)(Dr^] dv^^ 

(C-13) 

= Ju[(r+^U)'^(wJw‘) (r+*U) + 2(r+$U)'^WQW(*U). + 
(r+4U)'^WQW($u' )] dv + Juj^[(Dr^'^WQW' (Dr^ + 
2(Dr^’^wJw(Dr^ + ( Dr^'^WgW( 6 'r^] dv^ (C-14) 

The small motion equation is handled in a similar 
fashion, however it is complicated slightly by the 
deformation transformation matrix D. The product of the 
matrix terms involving D will result in scaler values, 
therefore the partial derivatives of D with respect to each 
small motion generalized coordinates will contribute 
separately. This requires the introduction of the operator g 
which represents a column vector composed of the matrix 
products, where the derivatives have been taken with respect 
to each small motion generalized coordinate. 

The symbol DK represents the deformation transformation 
matrix for which the partial derivative of matrix D have been 



d raKEi 

dtL a-J 
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taken with respect 


to 


each 


small 


motion generalized 


coordinate represented 


by 


K. 


This 


is demonstrated in 



equation (C-15) where K is 1. Applying the operator 9 we can 
express the scaler terms for the derivatives of D for each K, 
1-6 by equation (C-16). 




0 

0 

0 

0 



0 

0 

0 

0 



0 

0 

0 

0 



T T 

9 d"w^W Dj^ 



I- T T 
T T 

° “e"°K=2 

■ 

■ 

■ 

T T 



(C-15) 



(C-16) 



Separating the contributions from the link, and the load 
we present: 

For the link; 



9KE 

SfiT 



1 






4^w'^W(r+-^U)] dv 



(C-17) 



af[|rT’] = Ju[**^(w^w + w"^w)$u + + <t'^w^w$u + 

<t^(W^ U + w'^w“ ) (r+*U)] dv 



= Ju[-t^w'^W(r + *U) + $^w'^w<tu] dv 



(C-18) 

(C-19) 
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+ *'^w'^W$U + ^'^w'^ft*(c.+*U)] dv 



For the load: 



(C-20) 



?KE 






r " cJ'>’'S’'wD,(tJ dv^_ (C-21) 



K 



3KE^ 






'K 



tt**t •T« t«T«T •• 

r^D^(W W + W^W)Dj^r^+ r ^D^W^WDj^r^] dv^ 



T ■ T* T 



(C-22) 



9KE 






(C-23) 



^ ^ . !!!l ^ f r V„T {;5 ,V„T„ 5 . cVSWr,^ 

dtL aft J ar“i J^L»--L K ~L ~L K “L K-"L 

3U„ 3U„ 



t^D’'a-'«D„Ci,] dv,_ 



(C-24) 



Completing our development, the expressions for the 
potential energy are presented organized as the previous 
material and followed by the Lagrange equations. 



PE = j J u'^r'^CrU dx + jH(r + *U)'^w\ dv - J; 



dvj_^ 



(C-25) 



3PE ^ 

ae 



= - Ju(£l + <fU)'^wJa dv - JuLElJ^D'^WQa dVj 



(C-26) 
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Recalling equations (C-2), (C-23), and (C-26) we can make 
the appropriate substitutions and present the Lagrange 
equation for the large motion generalized coordinates. 



F = Ju [(EL+*U) WqW' ) (r+*U) + 2 ( E.+*U ) ( *5 ) + 

(EL+*U)'^WqW(^u’ )] dv + Juj^[(Dt^'^WQW"(Dr^ + 
2(Dr^'^WQW(6r^ + ( Dr^'^W^WC 6'r^] dv^ + 

- Ju(E- + *U)^WQg. dv - (C-27) 



We can simplify this expression further by the 
substitution of the integral terms defined in TABLE II and 
the separation of the second time derivative of the 
coordinate transformation matrix W by its' force 
contributions, equation (C-28). The contributions to 
Coriolis and centripetal forces are represented by Wr, termed 
the residual accelerations. The contribution to the general 
forces is represented by We. 



ft’ = W + W^ 
r 0 



(C-28) 



+u''l22(»X* "X’" * 

+ u’^l22(wJw)U* + HgCo’^W^W^D) + Hg(D’^W^WQD) + 2Hg(D'^wJwD) 



+ h^cd'^w^wS') - H^wjg. - u'^H2W^a ~ H^o'^wja 



(C-29) 
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Returning to the development o£ the small motion Lagrange 
equation the potential energy contribution remains to be 
presented. Again as in the previous case the contributions 
of the point load and the link will be developed separately. 
For the link: 



= f r'^CrU dx - dv 



(C-30 ) 



d r^KE-i 

dtL "tJ 
au 





(C-31) 



+ l22(w'^Wr.^^ ^1^ " ' (C-32) 

For the load: 

3PE. « m m m 

-^--KWadv^_ (C-33) 



dr3KE. 3KE ^ 9PE . „ 



(C-34) 



mm** TT*« TT** 

0 = Hg(DjJw^WD) + Hg(DjJw^WD ) + Hg(D^W^W Dj^) + 

Hg(D^w"WDj^) - H^DjJw a (C-35) 

Collecting terms and arranging the coefficients the two 
Lagrange equations can be written as the equations of motion 
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for the large and small generalized coordinates. It is these 
equations which must be solved by the computer simulation 
code . 

For the large motion: 



r T 








h2'"I"e''' " 




«’'l22'"I"> * 


_ 


«’'*2i'"I''e> * 


■ ■ 

e + 


g HjCd'^WqWDj^) 










L 


_ »5 









r 






2Ij2<wJS)0 - 



fji jp jp T HP 

Hi«ea t + H^o^w^a 



"’’i22<'‘X>“ ■ 

T T ■ ■ 

2H.(D^wiWD) + 

D o 



For the small motion: 



(C-36) 







l22(»^W) + 


l22“'\>“ * 


■ • 

e • + 


5 Hg(D^W^WDj^) 


g H^(dJw’'w^D) 










2 l22<w'^W) + 




• 




2 g Hg(Dj^w'^WDj^) 




u 


+ 



I22(W^W^) + + 

5 H 5 <Dyw^D,) 





U 




. 



[ - l2^(W^W^) + H2«\ + ? ] 



(C-37) 
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APPENDIX D 



LISTING OF COMPUTER SIMULATION CODE 



• THIS IS THE CONTROLLING PROGRAM FOR BOOM. IT DEFINES THE • 

• VARIABLES, DETERMINES THE COEFFICIENTS FOR THE EXPLICIT * 

• INTEGRATION, AND CALLS THE TWO SEPARATE INTEGRATOR OF THE 

« LARGE AND SMALL MOTION EQUATIONS. * 

• * 

• DEFINITION OF VARIABLES * 



A(2) H-D JOINT VARIABLE 

ALPHA(2) H-D JOINT VARIABLE * 

AREA CROSS SECTIONAL AREA OF BOOM 

E YOUNG'S MODULES 

EPSLON PARAMETER FOR INTEGRATION COEFFICIENTS * 

FINTIM FINISH TIME OF INTEGRATION * 

G SHEAR MODULES * 

GRAV GRAVITATIONAL FORCE (GLOBAL COORDINATES) * 

H TIME STEP USED IN THE INTEGRATION 

IOTA PARAMETER FOR INTEGRATION COEFFTCIENTS » 

lYY MOMENT OF INERTIA ABOUT THE Y AXIS 

IZZ MOMENT OF INERTIA ABOUT THE Z AXIS 

MU(1) MASS DENSITY OF BOOM (KG/M*»3) * 

MU(2) POINT MASS AT BOOM TIP (KG) 

Q LARGE MOTION GENERALIZED COORDINATE 9 

QDOT .LARGE MOTION ROTATIONAL VELOCITY 

T INDEXING VARIABLE TO PERMIT TIMING OF PRINT * 

TIME SEC * 

Q2DOT LARGE MOTION ROTATIONAL ACCELERATION 

TORK LARGE MOTION FORCING FUNCTION . » 

U(l) SMALL MOTION AXIAL DEFLECTION 

U(2) SMALL MOTION TRANSVERSE DEFLECTION (Y) 

U(3) SMALL MOTION TRANSVERSE DEFLECTION (Z) 

U(4) SMALL MOTION TORSIONAL DEFLECTION 

U(5) SMALL MOTION ROTATION ABOUT LOCAL Y AXIS * 

U(6) SMALL MOTION ROTATION ABOUT LOCAL Z AXIS * 



INTEGER PT,T 

REAL* 8 A( 2 ) , D( 2 ) , THETA ( 2 ), ALPHA ( 2 ), AREA, E,G, 
+ lYY, IZZ, MU(2), GRAV, 

+ U(6),UDOT(6),U2DOT(6), 

+ - Q,QDOT,Q2DOT, 

+ EPSLON, IOTA, H, TIME, FINTIM, TORK 



COMMON /BLOCKl/ A, D, THETA, ALPHA, AREA, E,G, I YY, IZZ, MU, GRAV 

COMMON /BLK15/ Q,QDOT,Q2DOT 

COMMON /BLK16/ U,UDOT,U2DOT 

GRAV=0 . DO 

TORK=100.D0 

MU(2)=44. 1 

MU(2)=0.D0 

FINTIM =2. DO 

H = .0006D0 

PT = 34 

EPSLON = 0.25D0 
IOTA = 0 . 5D0 
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A(1)=0. 0 
D( 1)=0 . 

THETA ( 1 )=Q 

ALPHA( 1)=2. *ATAN( 1. ) 

A(2)=5.486 

D(2)=0. 

ALPHA(2)=0. 

AREA=1. 9E-3 
E=70.E9 
G=26 . 0E9 
I?y=34. 63E-6 
122=34. 63E-6 
MU( 1)=2710.D0 



Q = O.DO 
QDOT = O.DO 
Q2DOT = O.DO 
DO 5 1=1,6 
U(I)=O.DO 
UDOT( I ) =0 . DO 
5 U2D0T( I ) =0 . DO 
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THETA ( 2 ) =0 . DO 
T=0 



DO 10 TIME=0,FINTIH.H 
T=T<-1 



CALL CCOEF(EPSLON,IOTA,H) 

CALL INTSHL(TIME,TORK) 

CALL INTLRG(TIME,H,TORK) 

IF( TIME. EQ. 0 . DO )CALL RINTKTIMEJ 
IF(T.EQ.PT)THEN 



CALL RINTKTIHE) 
T=0 



ENDIF 

CONTINUE 



END 

SUBROUTINE INTSHL( TIME, TORK) 



* THIS SUBROUTINE IMPLICITLY INTEGRATES THE SMALL MOTION EQUATION * 

* ♦ 

* ADDITIONAL DEFINITION OF VARIABLES * 



* Al-7 THE COEFFICIENTS FOR NUMERICAL INTEGRATION 

* FN THE SMALL MOTION FORCE VECTOR 

* FNN REFORMED SMALL MOTION FORCE VECTOR * 

* FQ THE LARGE MOTION FORCE VECTOR 

* GN GYROSCOPIC MATRIX 

* KN STIFFNESS MATRIX 

MN REFORMED SMALL MOTION INERTIA AND COUPLING * 

* MATRIX ♦ 

* MNN SMALL MOTION INERTIA MATRIX 

* MQN COUPLING INERTA MATRIX * 

* MQQ LARGE MOTION INERTIA MATRIX 

* TEMPl-5 TEMPORATY MATRICES FOR INTERMEDIATE VALUES • 

* PRD DUMMY ARGUMENT USED IN TROUBLE SHOOTING 

* WKAREA STORAGE FOR IMSL LEQT2F ROUTINE 

* (USED FOR LINEAR EQUATION SOLVING) 
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REAL*8 MQQ,MQN(6),MNN(6,6),GN(6, 6) ,KN(6.6),FN(6),FQ. 

+ U(6),UDOT(6),U2DOT(6), 

-J* MN(6^6) ,FNN(6) , 

+ A0.A1.A2,A3, A4, A5,A6, A7, 

+ WKAREA(54) ,PRD(6,6) ,TIME,TORK, 

+ TEMPI ( 6 ) , TEMP2 ( 6 ) , TEMP3 ( 6 ) , TEMP4 ( 6 ) , TEMPS ( 6.6) 

COMMON /BLK16/ U.UDOT.U2DOT 

COMMON /BLOCKS/ MQQ. MQN . MNN . GN . KN. FN. FQ 

COMMON /BLOCK?/ MN.FNN 

COMMON /BLOCKS/ AO . Al, A2. A3. A4. AS. A6. A? 

CALL TRANSF(TIME) 

CALL WMATRX 
CALL WDERIV 

IF ( TIME . EQ . 0 . DO ) CALL LRGCOF ( TORK ) 

CALL SMLCOF 
CALL REFORM^ 

DO 10 1=1.6 

TEMPKD = A0*U(I) + A2*UDOT(I) + A3*U2DOT(I) 

10 TEMP2(I) = A1*U(I) + A4*UDOT(I) + AS-'U2DOT(I) 

CALL MATMUL(PRD.MN.6.6.TEMP1.6. 1.TEMP3) 

CALL MATMUL(PRD.GN.6.6.TEMP2.6. 1.TEMP4) 

DO 20 1=1.6 

TEMPKD = FNN(I) + TEMP3(I) + TEMP4 ( I ) 

TEMP2(I) = U(I) 

DO 20 J=1.6 

20 TEMPS(I.J) = KN(I.J) + A0*MN(I.J) + A1*GN(I.J) 

M = 1 
N = 6 
IDGT = 1 

CALL LEQT2F( TEMPS. M.N.N. TEMPI. IDGT. WK AREA. lER) 

DO 30 1=1.6 
30 U(I) = TEMPKD 

DO 40 1=1.6 

TEMP3(I) = U2DOT(I) 

U2DOT(I) = A0*(U(I) - TEMP2(D) - A2*UDOT(I) - A3*TEMP3(I) 
40 UDOT(I) = UDOT(I) + A6*TEMP3(I) + A7-'U2DOT(I) 



END 
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SUBROUTINE INTLRG ( TIME. H, TORK) 



THIS SUBROUTINE INTEGRATES THE LARGE MOTION EQUATION EXPLICITLY * 



* ADDITIONAL DEFINITION OF VARIABLES * 

* * 



* 


Al-7 




INTEGRATION 


* 


* 


FN 


THE SMALL MOTION FORCE VECTOR 




* 




FQ 


THE LARGE MOTION FORCE VECTOR 




it 


* 


FQN 


REFORMED LARGE MOTION FORCE 




* 


It 


GN 






it 


* 


KN 


STIFFNESS MATRIX 




it 


* 


MN 


REFORMED SMALL MOTION INERTIA , 


AND COUPLING 


it 


* 




MATRIX 




* 


★ 


MNN 


SMALL MOTION INERTIA MATRIX 




* 


* 


MQN 


COUPLING INERTA MATRIX 




* 


* 


MQQ 


LARGE MOTION INERTIA MATRIX 




* 


it 


TEMPI 


TEMPORATY MATRIX FOR INTERMEDIATE VALUES 


* 


* 


PRD 


DUMMY ARGUMENT USED IN TROUBLE 


SHOOTING 


* 



* 



************************************************************************ 



REAL * 8 MQQ , MQN ( 6 ) , MNN (6,6),GN(6,6),KN(6,6),FN(6),FQ, 
+ U(6),UDOT(6),U2DOT(6), 

+ Q,QDOT,Q2DOT, 

+ A0,A1,A2,A3,A4,A5,A6,A7, 

+ H, TEMPI, PRD(6,l),FQN,TORK, TIME 

COMMON /BLOCKS/ MQQ, MQN, MNN, GN, KN, FN, FQ 
COMMON /BLK15/ Q,QDOT,Q2DOT 
COMMON /BLK16/ U,UDOT,U2DOT 
COMMON /BLOCKS/ AO , A1 , A2, A3 , A4, A5 , A6 , A7 

Q = Q + H’QDOT + A3*Q2DOT/AO 
QDOT = QDOT + A6»Q2DOT 

CALL TRANSF(TIME) 

CALL WMATRX 
CALL WDERIV 
CALL LRGCOF(TORK) 

TEMPI =0.00 
DO 10 1=1,6 

10 TEMPI = TEMPI + MQN(I)»U2DOT(I) 

FQN = FQ - TEMPI 

Q2DOT = (l.DO/MQQ)*FQN 
QDOT = QDOT + Q2DOT*A7 
Q = Q + Q2DOT/AO 
END 
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SUBROUTINE LRGCOF(TORK) 

★♦♦A******************************************************************** 



* THIS SUBROUTINE CALCULATES VALUES FOR COEFFICIENT MATRICES FOR * 

* LINK ONE. * 

* * 

* ADDITIONAL DEFINITION OF VARIABLES * 



Al-7. . 
CONST. 

DC 

DCT. . . 
DDOT. . 
DDOTT. 
DFORM. 
DFORMT 

FN 

FQ 

FQN. . . 

GN 

KN 

MAT*** 



MN. . . 

MNN. . 
MQN. . 
MQQ. . 
TEMPI 
ONE. . 



PRD. . . 

W 

WPRM. . 
WDOT. . 
WRESID 
Wl-14. 

XL 

XU 

Y.Z. . . 



THE COEFFICIENTS FOR NUMERICAL INTEGRATION * 

SINGLE VALUED TERM USED IN SUMMATION FOR 
FQ VALUE * 

DERIVATIVE OF DFORM MATRIX WRT SMALL MOTION * 

COORDINATE OF INTEREST * 

TRANSPOSE OF DC MATRIX * 

THE RATE OF DEFORMATION TRANSFORMATION MATRIX* 
THE TRANSPOSE OF THE DDOT MATRIX 
THE DEFORMATION TRANSFORMATION MATRIX * 

THE TRANSPOSE OF THE DFORM * 

THE SMALL MOTION FORCE VECTOR * 

THE LARGE MOTION FORCE VECTOR * 

REFORMED LARGE MOTION FORCE 

GYROSCOPIC MATRIX * 

STIFFNESS MATRIX * 

VARIOUS DIMENSIONED MATRICES OF INTERMEDIATE* 
VALUES {***) FOR THE VARIOUS SIZES * 

REFORMED SMALL MOTION INERTIA AND COUPLING * 

MATRIX * 

SMALL MOTION INERTIA MATRIX * 

COUPLING INERTA MATRIX * 

LARGE MOTION INERTIA MATRIX * 

TEMPORATY MATRIX FOR INTERMEDIATE VALUES ** 

ARRAY OF UNITY NECESSARY IN OBTAINING 
DERIVATIVE OF DFRM MATRIX . * 

DUMMY ARGUMENT USED IN TROUBLE SHOOTING * 

HARTENBERG DENAVIT TRANSFORMATION MATRIX * 

DERIVATIVE OF W WRT 9 * 

DERIVATIVE OF W WRT TIME * 

RESIDUAL ACCELERATION TRANSFORMATION * 

VARIOUS KERNELS FOR THE BOOM VOLUME INTEGRAL * 
LOWER LIMIT OF BOOM VOLUME INTEGRATION * 

UPPER LIMIT OF BOOM VOLUME INTEGRATION * 

DUMMY ARGUMENTS * 



* 

A***************************************************** 



INTEGER INDEX, NR, NC,NR2,NC2, TERM 

REAL*8 A(2) ,D( 2) ,THETA(2) , ALPHA(2) , AREA,E,G, 

+ IYY,IZZ,MU(2),GRAV, 

+ U ( 6 ) , UDOT ( 6 ) , U2DOT( 6 ) , 

+ UT(6), 

+ Q,QDOT,Q2DOT, 

+ W(4,4), 

+ WPRM ( 4, 4 ) ,WDOT( 4,4) , WRESID( 4, 4) , 

+ WD1(4,4),WD2(4,4),WD3(4,4), WD4 (4,4), WD6 (4,4), 

+ WD9(4,4) ,WD10(4,4) ,WD13(4,4) ,WD14(4,4) , 

+ MQQ,MQN(6) ,MNN(6,6),GN(6,6) ,KN(6, 6) ,FN(6) ,FQ, 

+ MAT(6) ,MAT4A(4) ,MAT4B(4) ,MAT64A(6,4) ,MAT64B(6, 4) ,MATPRD(6,6) , 

+ CONST,Z(6,6), Y(6,6),XL,XU,ONE(6) ,MAT44A(4,4),MAT44B(4,4) , 

+ PTMASS(4) ,GRAVEC(4) ,WT(4,4) ,WPRMT(4,4) ,PRD(6,6) ,DC(4, 4) , 

+ DFRM(4,4) ,DFRMT(4,4) ,DDOT(4,4) ,DDOTT( 4, 4) ,DCT(4,4) ,TORK 
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COMMON /BLOCKl/ A, D, THETA. ALPHA , AREA. E, G, I YY, IZ2, MU . GRAV 
COMMON /BLK15/ Q.QDOT.Q2DOT 
COMMON /BLK16/ U.UDOT.U2DOT 
COMMON /BLOCK3/ W 

COMMON /BLOCK4/ WDl. WD2 . WD3. WD4 . WD6. WD9. WDIO . WD13. WD14 
COMMON /BLOCKS/ MQQ. MQN. MNN. GN. KN. FN. FQ 
COMMON /BLOCKS/ WPRM. WDOT. WRESID 

XL=-A(2) 

XU=0.D0 

MQQ =0.00 
FQ=T0RK 
DO 3 1=1.6 

3 MQN(I)=0.D0 

DO 5 J=1.4 

GRAVEC( J)=0.D0 
5 PTMASS(J)=0.D0 

GRAVEC(4)=GRAV 

PTMASS(1)=MU(2) 

TERM=0 

DO 7 J=1.6 

7 ONE(J)=1.DO 

CALL DFORM( 1.7.U.DFRM) 

CALL DFORM(0. 7.UDOT.DDOT) 

CALL TRANSP(TERM.DFRM.4.4.DFRMT) 

CALL TRANSP(TERM.DDOT.4.4.DDOTT) 



CALL TRANSP(TERM.U.6. l.UT) 

CALL TRANSP(TERM.W.4.4.WT) 

CALL TRANSP(TERM.WPRM.4.4.WPRMT) 



*************************inf»**im*'0n*itimitimitinti,i,niti,ttitit**mt*innm**iiiiiri,t,i, 

CREATING THE COEFFICIENT MATRIX MQQ 



INDEX = 11 

CALL BINT ( TERM . XL. XU . Z . Y. 1. 1 . INDEX. WD6 . CONST ) 
HQQ=MQQ+CONST 

INDEX = 12 

CALL BINT(TERM.XL.XU.Z. Y. 1.6. INDEX. WD6 . HAT ) 
CALL MATHUL( PRD. HAT. 1. 6. U. 6 . 1 . CONST ) 
MQQ=MQQ+CONST 

INDEX = 21 

CALL BINT(TERM.XL.XU.Z.Y.6.1.INDEX.WD6.HAT) 
CALL MATHUL(PRD.UT. 1.6.HAT.6. l.CONST) 
MQQ=MQQ+CONST 

INDEX = 22 

CALL BINT(TERM.XL.XU.Z. Y.6.6.INDEX.WD6.MATPRD) 
CALL MATMUL(PRD.MATPRD,6.6.U.6. l.MAT) 

CALL MATMUL( PRD. UT. 1. 6. MAT. 6. l.CONST) 
HQQ=HQQ+CONST 

THE CONTRIBUTION OF THE POINT LOAD 

CALL HATMUL( PRD . WD6 .4.4. DFRM .4.4. H AT4 4 A ) 

CALL HATHUL(PRD.DFRMT.4.4.MAT44A.4. 4.HAT44B) 



ONLY ELEMENT MAT44B(1.1) WILL SURVIVE FOLLOWING POST MULTIPLICATION BY 
LOCAL POSITION VECTOR AND PRE MULTIPLICATION BY ITS TRANSPOSE. 



C0NST=MU(2) •HAT44Ba.l) 

MQQ=MQQ+CONST 

n*it*imit*it************************************************************** 

* CREATING THE COEFFICIENT MATRIX MQN 

ifk***********************************************'^********************* 

INDEX=12 

CALL BINT(TERM,XL,XU,Z,Y, 1,6, INDEX, WD2, MAT) 

DO 8 1=1,6 

8 MQN(I)=MQN(I) + MAT(I) 

INDEX = 22 

CALL BINT(TERM,XL,XU, Z,Y,6,6,INDEX, WD2,MATPRD) 

CALL MATMUL(PRD,UT, 1,6,MATPRD,6,6,MAT)- 
DO 9 1=1,6 

9 MQN(I) = MQN(I) + MAT(I) 

* THE CONTRIBUTION OF THE POINT LOAD 



DO 10 K=l,6 

CALL DFORM(0,K,ONE,DC) 

CALL MATMUL(PRD,WD2,4, 4,DC,4,4,MAT44A) 

CALL MATMUL(PRD, DFRMT, 4 , 4 , MAT44 A, 4 , 4, MAT4 4B) 

* ONLY ELEMENT MAT44B(1,1) WILL SURVIVE FOLLOWING POST MULTIPLICATION BY 

* LOCAL POSITION VECTOR AND PRE MULTIPLICATION BY ITS TRANSPOSE. 

10 MQN(X)=MQN(K) + MU ( 2 ) *MAT44B ( 1 , 1 ) 



*********************************************************************** 

• CREATING THE COEFFICIENT MATRIX FQ 



INDEX = 11 

CALL BINT(TERM,XL,XU,Z, Y, 1, 1,INDEX, WD14,CONST) 
FQ=FQ -CONST 

INDEX = 12 

CALL BINT (TERM, XL, XU, Z,Y, 1,6, INDEX, WDIO, MAT) 
CALL MATMUL(PRD,MAT, l,6,UDOT,6, 1, CONST) 

FQ=FQ-2. ‘CONST 

CALL BINT(TERM,XL,XU,Z,Y, 1,6,INDEX, WD14,MAT) 
CALL MATMUL(PRD,MAT, 1, 6 , U, 6, 1 , CONST) 

FQ=FQ-CONST 

INDEX = 21 

CALL BINT(TERM,XL,XU,Z, Y, 6, 1, INDEX, WD14, MAT ) 
CALL MATMUL(PRD,UT, 1,6,MAT,6, 1, CONST) 
FQ=FQ-CONST 

INDEX = 22 

CALL BINT(TERM,XL,XU,Z, Y, 6, 6 , INDEX, WDIO , MATPRD) 
CALL MATMUL(PRD,MATPRD,6,6,UDOT,6, 1,MAT) 

CALL MATMUL(PRD,UT, 1,6,MAT,6, l,CONST) 

FQ=FQ-2. ‘CONST 

CALL BINT ( TERM, XL, XU, Z, Y, 6, 6, INDEX, WDl 4, MATPRD) 
CALL MATMUL(PRD,MATPRD,6,6,U,6, 1,MAT) 

CALL MATMUL(PRD,UT, 1,6,MAT,6, 1, CONST) 
FQ=FQ-CONST 

‘ THE CONTRIBUTION OF THE POINT LOAD 
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CALL HATMULCPRD, WD14,4.4,DFRH.4,4.MAT44A) 

CALL MATMUL( PRD. DFRMT, 4 . 4. MAT44A. 4 . 4. MAT44B ) 

* ONLY ELEMENT MAT44B(1.1) WILL SURVIVE FOLLOWING POST MULTIPLICATION BY 

* LOCAL POSITION VECTOR AND PRE MULTIPLICATION BY ITS TRANSPOSE. 

CONST=MU(2) *MAT44B( 1. 1) 

FQ=FQ-CONST 

CALL MATMULCPRD. WDIO. 4.4.DDOT. 4. 4.MAT44A) 

CALL MATMULC PRD.DFRMT, 4. 4.MAT44A. 4. 4.MAT44B ) 

* ONLY ELEMENT HAT44B(1.1) WILL SURVIVE FOLLOWING POST MULTIPLICATION BY 

* LOCAL POSITION VECTOR AND PRE MULTIPLICATION BY ITS TRANSPOSE. 

CONST=2 . D0*MU ( 2 ) *MAT44B( 1. 1 ) 

FQ=FQ-CONST 

* THE CONTRIBUTION OF POTENTIAL ENERGY 

INDEX = 10 

* HI. WDl WILL BE CARRIED AS A DUMMY ARGUMENT 

CALL BINTCTERM.XL.XU.2. Y. 1. 4. INDEX. WDl. MAT4A ) 

CALL MATMULCPRD.MAT4A. 1. 4. WPRMT. 4. 4 . MAT4B ) 

CALL MATMULCPRD.MAT4B. 1. 4.GRAVEC. 4. 1. CONST) 

FQ=FQ+CONST 

INDEX = 20 

* H2. WDl WILL BE CARRIED AS A DUMMY ARGUMENT 

CALL BINT ( TERM. XL. XU. 2. Y. 6. 4. INDEX. WDl. MAT64 A) 

CALL MATMUL(PRD.UT.1.6.MAT64A.6.4.MAT4A) 

CALL MATMUL(PRD.MAT4A.1.4.WPRMT.4.4.MAT4B) 

'CALL MATMUL(PRD.MAT4B.1.4.GRAVEC.4.1.CONST) 

FQ=FQ+CONST 

* THE CONTRIBUTION OF THE POINT LOAD TO THE POTENTI/^L ENERGY 

CALL MATMUL( PRD. DFRMT. 4. 4. WPRMT. 4. 4. MAT44A ) 

CALL MATMUL(PRD.MAT44A.4.4.GRAVEC.4/1.MAT4B) 

CALL MATMUL(PRD.PTMASS. 1.4.MAT4B. 4. 1. CONST) 

FQ=FQ+CONST 



RETURN 

END 
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SUBROUTINE SMLCOF 

************************************************************************ 



* THIS SUBROUTINE CALCULATES VALUES FOR COEFFICIENT MATRICES FOR * 

* LINK ONE. * 

* * 

* ADDITIONAL DEFINITION OF VARIABLES * 

* * 



* 



* 



Al-7. . 
CONST. 

DC 

DCT. . . 
DDOT. . 
DDOTT. 
DFORM. 
DFORHT 

FN 

FQ 

FQN. . . 

GN 

KN 

HAT*** 



HN 

MNN. . . 
HQN. . . 
HQQ. . . 
TEMPI. 
ONE. . . 

PRD. . . 

U 

WPRM. . 
WDOT. . 
WRESID 
Wl-14. 

XL 

XU 

Y,Z. . . 



THE COEFFICIENTS FOR NUMERICAL INTEGRATION * 

SINGLE VALUED TERM USED IN SUMMATION FOR * 

FQ VALUE * 

DERIVATIVE OF DFORM MATRIX WRT SMALL MOTION * 
COORDINATE OF INTEREST * 

TRANSPOSE OF DC MATRIX * 

THE RATE OF DEFORMATION TRANSFORMATION MATRIX* 
THE TRANSPOSE OF THE DDOT MATRIX * 

THE DEFORMATION TRANSFORMATION MATRIX * 

THE TRANSPOSE OF THE DFORM * 

THE SHALL MOTION FORCE VECTOR * 

THE LARGE MOTION FORCE VECTOR * 

REFORMED LARGE MOTION FORCE * 

GYROSCOPIC MATRIX * 

STIFFNESS MATRIX * 

VARIOUS DIMENSIONED MATRICES OF INTERMEDIATE * 
VALUES (***) FOR THE VARIOUS SIZES * 

REFORMED SMALL MOTION INERTIA AND COUPLING * 

MATRIX * 

SHALL MOTION INERTIA MATRIX 

COUPLING INERTA MATRIX * 

LARGE MOTION INERTIA MATRIX * 

TEMPORATY MATRIX FOR INTERMEDIATE VALUES * 

ARRAY OF UNITY NECESSARY IN OBTAINING 
DERIVATIVE OF DFRM MATRIX * 

DUMMY ARGUMENT USED IN TROUBLE SHOOTING * 

HARTENBERG DENAVIT TRANSFORMATION MATRIX * 

DERIVATIVE OF W WRT 9 * 

DERIVATIVE OF W WRT TIME * 

RESIDUAL ACCELERATION TRANSFORMATION 
VARIOUS KERNELS FOR THE BOOH VOLUME INTEGRAL * 
LOWER LIMIT OF BOOH VOLUME INTEGRATION * 

UPPER LIMIT OF BOOH VOLUME INTEGRATION * 

DUMMY ARGUMENTS * 



************************************************************************ 



INTEGER INDEX. NR. NC,NR2,NC2. TERM 

REAL*8 A ( 2 ) . D ( 2 ) . THETA ( 2 ) . ALPHA ( 2 ) . AREA. E. G, 

^ IYY.IZZ.MU(2).GRAV. 

-f W ( 4 4 ) 

+ U(6).UD0T(6).U2D0T(6). 

+ UT(6). 

+ Q.QDOT.Q2DOT. 

+ WPRM (4.4). WDOT (4.4). WRESID (4.4). 

WD1(4.4).WD2(4.4).WD3(4.4).WD4(4.4).WD6(4.4). 

+ WD9(4.4).WD10(4.4).WD13(4.4).WD14(4.4). 

+ MQQ.MQN(6).MNN(6.6).GN(6.6).KN(6.6).FN(6).FQ. 

+ HAT(6).HAT4A(4).HAT4B(4).HAT64A(6.4).HAT64B(6.4).HATPRD(6.6). 

CONST.Z(6.6). Y(6.6).XL.XU.ONE(6).MAT44A(4.4).HAT44B(4.4). 

+ PTMASS(4).GRAVEC(4).WT(4.4).WPRHT(4,4).PRD(6.6),DC(4.4). 

DFRH(4.4).DFRMT(4.4).DDOT(4.4).DDOTT(4,4).DCT(4.4) 

COMMON /BLOCKl/ A . D. THETA. ALPHA . AREA . E. G. I YY . IZZ . HU. GRAV 
COMMON /BLK15/ Q.QDOT.Q2DOT 
COMMON /BLK16/ U.UDOT.U2DOT 
COMMON /BLOCK3/ W 

COMMON /BLOCK4/ WDl . WD2. WD3 . WD4 . WD6 . WD9 . WDIO . WD13 . WD14 
COMMON /BLOCKS/ MQQ. HQN. MNN. GN. KN. FN. FQ 
COMMON /BLOCKS/ WPRM. WDOT. WRESID 
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3 



XL=-A(2) 

XU=O.DO 
DO 3 I«l,6 
FN(I)=0.D0 
DO 3 J*l,6 

GN(I, J)=0.D0 
KN(I,J)=»0.D0 
MNN(I,J)=0.D0 
DO 5 J=«l,4 

GRAVEC(J)=0.D0 
5 PTMASS(J)=0.D0 

GRAVED (4 )=GRAV 
PTMASS(1)=MU(2) 

TERM«0 
DO 7 J»l,6 
7 ONE(J)*1.DO 

CALL DFORM( 1,7,U,DFRM) 

CALL DFORM(0,7,UDOT,DDOT) 

CALL TRANSP(TERM,DFRM,4M,DFRMT) 
CALL TRANSP(TERM,DDOT,4,4,DDOTT) • 
CALL TRANSP(TERM,U,6, 1,UT) 

CALL TRANSP(TERM,W,4,4,MT) 

CALL TRANSP(TERM,WPRM,4,4,WPRMT) 

* CREATING THE COEFFICIENT MATRIX MNN 



INDEX =« 22 

CALL BINT(TERM,XL,XU,2,y,6,6,INDEX^«Dl,MATPRD) 

DO 12 1=1,6 
DO 12 J=l,6 

12 MNN(I,J)=MNN(I,J) + MATPRD(I,J) 

* THE CONTRIBUTION OF THE POINT LOAD 

DO 15 K=i,6 

CALL DFORM(0,X,ONE,DC) 

CALL TRANSP(TERM,DC,4,4,DCT) 

CALL MATMUL(PRD,DCT,4,4,MD1,4,4,MAT44A) 

CALL DFORM(0,K,ONE,DC) 

CALL MATMUL(PRD,MAT44A,4,4,DC,4M,MAT44B) 

* ONLY ELEMENT MAT44B(1,1) MILL SURVIVE FOLLOWING POST MULTIPLICATION BY 

* LOCAL POSITION VECTOR AND PRE MULTIPLICATION BY ITS TRANSPOSE. 



15 MNN(K,K)*MNN(K,K) ^ MU ( 2 ) *MAT44B( 1, 1) 

A********************************************************************** 

* CREATING THE COEFFICIENT MATRIX GN (UDOT) 

A********************************************************************** 



INDEX = 22 

CALL BINT(TERM,XL,XU,2. Y,6,6,INDEX,WD9,MATPRD) 

DO 22 1=1.6 
DO 22 J=l,6 

22 GN(I, J)=GN( I,J)+2*MATPRD( I.J) 

* THE CONTRIBUTION FROM THE POINT LOAD 

DO 30 J=1.6 

CALL DFORMCO. J.ONE.DC) 

CALL MATMUL(PRD.WD9.4.4.DC.4.4.MAT44A) 

CALL TRANSP(TERM.DC.4,4,DCT) 

CALL MATMUL( PRD.DCT. 4. 4.MAT44A, 4. 4.MAT44B ) 

* ONLY ELEMENT MAT44B(1.1) MILL SURVIVE FOLLOWING POST MULTIPLICATION BY 

* LOCAL POSITION VECTOR AND PRE MULTIPLICATION BY ITS TRANSPOSE. 
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30 



CONST=2»MU( 2)*MAT44B(1, 1) 
GN( J,J)=GN( J. J)+CONST 



• CREATING THE COEFFICIENT HATRIX KN(0) 

»*****««**»»«•*»»*•»»•***•*•**•*•*»•*•••••••••*••*••••••*•••••*•••*•*•• 



INDEX * 22 

CALL BINT(TERH,XL,XU,Z, Y,6,6,INDEX,WD13,MATPRD) 
DO 24 1*1,6 
DO 24 J*l,6 

24 KN( I,J)*KN(I,J)+MATPRD(I,J) 



* THE CONTRIBUTION OF POTENTIAL ENERGY FROM THE LINK 

* Kl, WDl WILL BE CARRIED AS A DUMMY ARGUMENT 

INDEX * 1 

CALL BINT(TERM,XL,XU,Z,Y, 6 , 6 , INDEX, WD1,MATPRD) 
DO 28 1-1,6 
DO 28 J=l,6 

28 KN( I,J)*KN(I,J)+MATPRD(I,J) 



* THE CONTRIBUTION FROM THE' POINT LOAD 



DO 32 J*l,6 

CALL DFORMCO, J,ONE,DC) • 

CALL MATMUL(PRD,WD13,4,4,DC,4,4,MAT44A) 

CALL DFORM(0, J,ONB,DC) 

CALL TRANSP(TERM,DC,4,4,DCT) 

CALL MATMUL(PRD,DCT,4,4,MAT44A,4,4,MAT44B) 

• ONLY ELEMENT MAT44B(1,1) WILL SURVIVE FOLLOWING POST MULTIPLICATION BY 

* LOCAL POSITION VECTOR AND PRE MULTIPLICATION BY ITS TRANSPOSE. 



CONST-MU(2)*MAT44B( 1,1) 

32 KN( J,J)=KN(J, J)+CONST 

A********************************************************************** 

• CREATING THE COEFFICIENT MATRIX FN 



INDEX * 21 

CALL BINT(TERM,XL,XU,Z, Y,6,1,INDEX,WD13,MAT) 

DO 20 J*l,6 

20 FN(J)*FN(J)-MAT(J) 

• THE CONTRIBUTION OF POTENTIAL ENERGY FROM THE LINK 

• H2, WDl WILL BE CARRIED AS A DUMMY ARGUMENT 

INDEX * 20 

CALL BINT(TERM,XL,XU,Z, Y,6,4,INDEX,WD1,MAT64A) 

CALL MATMUL(PRD,MAT64A, 6,4,WT,4, 4,MAT64B) 

CALL MATMUL(PRD,MAT64B,6,4,GRAVEC,4, 1,MAT) 

DO 26 J=l,6 

26 FN( J)=FN( J)+MAT( J) 

• THE CONTRIBUTION OF THE POINT LOAD TO THE POTENTIAL ENERGY 

DO 36 J=l,6 

CALL DFORM(0, J,ONE,DC) 

CALL TRANSP(TERM,DC,4,4,DCT) 

CALL MATMULC PRD, DCT, 4, 4 , WT, 4, 4 , MAT44A ) 

CALL MATMUL(PRD,MAT44A,4,4,GRAVEC,4, 1,MAT4A) 

CALL MATMUL(PRD,PTMASS,1,4,MAT4A,4, l,CONST) 

36 FN( J)=FN(J)+CONST 

RETURN 

END 
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SUBROUTINE WMATRX 

A********************************************************************** 



THIS SUBROUTINE MULTIPLIES THE HARTENBERG-DENAVIT MATRICES * 

TOGETHER TO OBTAIN THE COORDINATE TRANSFORMATION MATRIX * 

* 

ADDITIONAL DEFINITION OF VARIABLES 



* 

HD THE HARTENBERG DENAVIT TRANSFORMATION MATRIX * 

W TRANSFORMATION MATRIX * 



A********************************************************************** 



REALMS W(4,4),HD(2,4,4),WSUM(4,4),SUM 

COMMON /BLOCK2/ HD 
COMMON /BLOCK3/ M 

DO 10 L»l,4 

DO 10 M=l,4 

WSUM(L,M)=0.D0 

10 W(L,M)=0.D0 

DO 40 J=^l,4 
DO 30 K=l,4 
SUM=0.0 
DO 20 M«l,4 

SUM « SUM + HD(1,J,M) * HD(2,M,K) 

20 CONTINUE 

WSUM(J,K)*SUH 
30 CONTINUE 
40 CONTINUE 

DO 50 I»l,4 
DO 50 J=l,4 

50 W(I,J)=WSUM(I,J) 

END 

SUBROUTINE WDERIV 

A********************************************************************** 



THIS SUBROUTINE CREATES THE DERIVATIVES OF THE COORDINATE MATRIX * 
W WRT Q THE LARGE MOTION GENERALIZED COORDINATE 9 WPRM * 

AND TIME MDOT AND WRESID. * 

* 

ADDITIONAL DEFINITION OF VARIABLES * 



* 

HD THE HARTENBERG DENAVIT TRANSFORMATION MATRIX ♦ 

HDPRM DERIVATIVE OF HD( 1) MRT 9 * 

W TRANSFORMATION MATRIX • 

WT TRANSPOSE OF THE TRANSFORMATION MATRIX • 

WPRM DERIVATIVE OF W WRT 9 * 

WPRMT TRANSPOSE OF THE DERIVATIVE OF W WRT 9 * 

WDOT DERIVATIVE OF W WRT TIME * 

WDOTT TRANSPOSE OF THE DERIVATIVE OF W WRT TIME * 

WRESID RESIDUAL ACCELERATION TRANSFORMATION * 

WRSIDT TRANSPOSE OF THE RESIDUAL ACCELERATION * 

TRANSFORMATION 

Wl-14 VARIOUS KERNELS FOR THE BOOM VOLUME INTEGRAL * 



A********************************************************************** 



INTEGER TERM 

REALMS W(4.4),HD(2,4,4), WPRM (4,4), WDOT (4,4), WRESID ( 4,4), 

+ Q,QDOT,Q2DOT, 

+ WT ( 4 , 4 ) , WPRMT ( 4,4), WDOTT ( 4 , 4 ) , WRSIDT ( 4 , 4 ) , 

+ HDPRM (4,4),HD2(4,4), PRD (4,4), 

+ WD1(4,4),WD2(4,4),WD3(4,4),WD4(4,4),WD6(4,4),WD9(4,4), 

+ WD10(4,4) ,WD13(4,4),WD14(4,4) 
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COMHON /BLK15/ Q,QDOT,Q2DOT 
COMMON /BLOCK2/ HD 
COMMON /BLOCKS/ U 

COMMON /BLOCK4/ WDl, WD2,WD3, WD4, WD6, WD9,WD10, WD13, WD14 
COMMON /BLOCKS/ WPRM, WDOT, WRESID 

DO 5 K=l,4 
DO 5 L=l,4 

HDPRM(K,L)=0.0 

HD2(K,L)=HD(2,K.L) 

WDOT(K,L)=0.0 

5 WRESID(K,L)=0.0 

DO 10 K « 1,4 

HDPRM(2,K)=-HD(1,3,K) 

10 HDPRM(3,K)= HD(1,2,K) 

CALL M ATMUL ( PRD, HDPRM , 4 , 4 , HD 2 , 4 , 4 , WPRM ) 

* THE FOLLOWING MATRIX HAS A QDOT COEFFICIENT 

DO 15 1=1,4 

DO 15 J=l,4 

IF(ABS(WPRM(I,J) ) .LE.l.E-30) THEN 
WDOT(I,J)=O.DO 

ELSE 

WDOT(I,J)=WPRM(I, J)*QDOT 
ENDIF 

15 CONTINUE 

» CREATING THE RESIDUAL ACCELERATION TRANSFORMATION MATRIX 
* THIS MATRIX HAS A QDOT**2 COEFFICIENT 

DO 20 K * 1,4 

IF (ABS(WDOT( 2, K) ) .LE.l.E-30) THEN 
WRESID(2,K)=0.D0 

ELSE 

WRESID ( 2 , K ) =-WDOT ( 3 , K ) *QDOT 
ENDIF 

IF ( ABS ( WDOT ( 3 , K ) ) . LE . 1 . E-30 ) THEN 
WRESID(3,K)=O.DO 

ELSE 

WRESID ( 3 , K ) =WDOT ( 2 , K ) *QDOT 
ENDIF 

20 CONTINUE 

CALL TRANSP(TERM,W,4,4,WT) 

CALL TRANSP(TERM,WPRM,4,4,WPRHT) 

CALL TRANSP(TERM,WDOT,4,4,WDOTT) 

CALL TRANSP( TERM, WRESID, 4,4, WRSIDT) 

CALL MATMUL (PRD,WT,4,4,W,4,4,WD1) 

CALL MATMULC PRD , WPRMT, 4,4,W,4,4,WD2) 

CALL MATMUL( PRD, WDOTT, 4 , 4 , W, 4 , 4 , WD3 ) 

CALL HATHUL( PRD, WRSIDT, 4, 4 , W, 4 , 4 , WD4 ) 

CALL MATMUL(PRD, WPRMT,4,4,WPRM, 4,4,WD6) 

CALL MATMUL( PRD, WT, 4 , 4 , WDOT, 4,4, WD9 ) 

CALL MATMUL( PRD, WPRMT, 4,4, WDOT, 4,4, WDIO ) 

CALL HATMUL( PRD, WT, 4, 4 , WRESID, 4 , 4 , WD13 ) 

CALL MATMUL( PRD, WPRMT, 4 , 4, WRESID, 4 , 4, WD14 ) 

END 
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SUBROUTINE INTPRD ( TERM, INDEX, WD, X, ROW, COU PROD ) 



THIS SUBROUTINE CREATES THE PRODUCT ARRAY OF THE SHAPE FUNCTIONS • 
AND THE COORDINATE TRANSFORMATION MATRICES W(I) 



INTEGER INDEX, ROW, COL, TERM 

REALMS PROD(ROW,COL),R(4),P(4,6),WD(4,4), 

+ A ( 2 ) , D( 2 ) , THETA( 2 ) , ALPHA( 2 ) , AREA, E, G, 

+ IYY,IZZ,MU(2),GRAV, 

+ P1(4,6),P1T(6,4),GMAT(4,6),GT(6,4),G1(4,6),G1T(6,4), 

+ WMAT(4,4),C(4,4),RC(4,1), 

+ PRD(6,6),X 

COMMON /BLOCXl/ A, D, THETA, ALPHA, AREA, E, G, I YY, IZZ , MU, GRAV 

IF (INDEX .EQ. 11) THEN 
CALL VECTR(X,R) 

PROD(l,l)»WD(l,l) + R(2)’*WD(2,1) + R(2)*WD(1,2) 

PROD(l,l)=PROD(l,l)+R(2)**2*WD(2,2)+IYY*WD(3,3)/(MU(l)*AREA) 

PROD(l,l)*PROD(l,l)+IZZ*WD(4,4)/(MU(l)*AREA) 

END IF 

IF (INDEX .EQ. 12) THEN 
CALL VECTR(X,R) 

CALL C0EFP1(X,P1) 

PROD(l,l)=WD(l,2) •Pl(2,l)+R( 2)*WD(2,2) *P1(2,1) 

PROD( 1,2)=WD( 1,3)*P1(3,2)+R(2) *WD(2,3)*P1(3,2) 
PROD(l,3)=WD(l,4)*Pl(4,3)+R(2)*WD(2,4)*Pl(4,3) 
PROD(l,4)=(-I2Z*WD(4,3) + I YY*WD( 3, 4 ) ) / ( MU(.l ) *AREA ) 
PROD(l,5)=WD( 1,4) *P1(4,5)+R(2)*WD(2,4) *P1(4,5) 
PROD(l,6)=WD(l,3)*Pl(3,6)+R(2)*WD(2,3)*Pl(3,6) 

END IF 



IF (INDEX .EQ. 21) THEN 
CALL C0EFP1(X,P1) 

CALL VECTR(X,R) 

PR0D(1,1)=WD(2,1)»P1(2,1)+R(2)*WD(2,2)*P1(2,1) 
PR0D(2,1)=WD(3,1)*P1(3,2)+R(2)*WD( 3,2) *P1(3,2) 
PROD(3,l)=WD(4,l)*Pl(4,3)+R(2)*WD(4,2)*PT(4,3) 
PROD(4,l)-(-IZZ*WD(3,4) + I YY*WD ( 4 , 3 ) ) / ( MU ( 1 ) * AREA ) 
PROD(5,l)=WD(4,l) *P1(4,5)+R( 2)*WD(4,2) *P1(4,5) 
PROD(6,l)«WD(3, 1) *P1(3,6)+R(2)*WD(3,2)*P1(3,6) 

END IF 



IF (INDEX .EQ. 22) THEN 



CALL C0EFP1(X,P1) 
PROD(l,l)=Pl(2,l) * 
PR0D(1,2)»P1(2,1) * 
PR0D(1,3)*P1(2,1) * 
PROD(1,4)»O.DO 
PR0D(1,5)-P1(2,1) * 
PR0D(1,6)=P1(2,1) * 
PROD(2,l)=PROD(l,2) 
PR0D(2,2)=P1(3,2) * 
PR0D(2,3)=P1(3,2) * 
PROD(2,4)=O.DO 
PR0D(2,5)=P1(3,2) * 
PR0D(2,6)=»P1(3,2) * 
PROD(3,l)=PROD(l,3) 
PROD(3,2)=PROD(2,3) 
PR0D(3,3)=P1(4,3) * 
PROD(3,4)*O.DO 
PR0D(3,5)=P1(4,3) * 
PR0D(3,6)=P1(4,3) * 
PROD(4,1)=O.DO 
PROD(4,2)*O.DO 
PROD(4,3)=O.DO 



WD(2,2) 


* 


Pl(2,l) 


WD(2,3) 


* 


Pl(3,2) 


WD(2,4) 


* 


Pl(4,3) 


WD(2,4) 


* 


Pl(4,5) 


WD(2,3) 


* 


Pl(3,6) 


WD(3,3) 


* 


Pl(3,2) 


WD(3,4) 


* 


Pl(4,3) 


WD(3,4) 


* 


Pl(4,5) 


WD(3,3) 


* 


Pl(3,6) 


WD(4,4) 


* 


Pl(4,3) 


WD(4,4) 


* 


Pl(4,5) 


WD(4,3) 


* 


Pl(3,6) 
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PROD(4,4)»(IZZ • HD(3,3) + I?Y • «D( 4, 3 ) ) / (MU( 1 ) 'AREA ) 

PROD(4,5)-O.DO 

PROD(4,6)»O.DO 

PROD(5,l)=PROD(l,5) 

PROD(5,2)=PROD(2,5) 

PROD(5,3)-PROD(3,5) 

PROD(5,4)=O.DO 

PR0D(5,5)=.P1(4,5) • «D(4,4) * Pl(4,5) 

PR0D<5,6)=P1(4,5) • HD(4,3) • Pl(3,6) 

PR0D(6, 1)=PR0D( 1,6) 

PROD(6,2)»PROD(2,6) 

PR0D(6,3)»PR0D(3,6) 

PROD(6.4)=O.DO 

PROD(6,5)=PROD(5,6) 

PR0D(6,6)»P1(3,6) • «D(3,3) • PK3.6) 

END IF 

IF (INDEX ,EQ. 1) THEN 
CALL COEFGKX.Gl) 

CALL TRANSP(TERM,G1,4,6,G1T) 

CALL STIFNS(C) 

CALL MATMUL(PRD,C,4,4,G1,4,6,G1) 

CALL MATMUL(PRD,G1T,6,4,G1.4,6,PR0D) 

END IF 

IF (INDEX .EQ. 10) THEN 
CALL VECTR(X.R) 

CALL TRANSP(TERM,R, 4,1, PROD) 

END IF 

IF (INDEX .EQ. 20) THEN 
CALL C0EFP1(X,P1) 

CALL TRANSP(TERM,P1,4,6,PR0D) 

END IF 
END 



SUBROUTINE COEFGl (X,G1) 



* THIS SUBROUTINE CALCULATES THE VALUES FOR THE GAMAl COEFFICIENT • 

* • 

********************************<»*************************************** 



REAL*8 A(2),D(2),THETA(2),ALPHA(2),AREA,E,G, 

+ IYY,I2Z,MU(2),GRAV, 

+ G1(4,6),L,X,K 

COMMON /BLOCXl/ A, D, THETA, ALPHA, AREA, E,G, lYY, I2Z, MU, GRAV 

L=A(2) 

X = X/L 



DO 10 J=1^4 
DO 10 M=l,6 

10 G1(J,M)=0.D0 



Gl(l,l) 

Gl(2,3) 

Gl(2,5) 

Gl<3,2) 

Gl(3,6) 

Gl(4,4) 

END 



1/L 

-(12. *X + 6. )/L*»2. 
(6. »X + 4, )/L 
GK2.3) 

Gl(2,5) 

GKl, 1) 
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DO 10 


J = l,4 


S(I,J)-0.D0 


P = lYVtIZZ 


S(l,l) 


= E-AREA 


S(2,2) 


- E*IYY 


SO, 3) 


= E*IZZ 


S(4,4) 


= G*P 



SUBROUTINE STIFNS (S) 

* * 

* THIS SUBROUTINE CREATES THE STIFFNESS MATRICES FOR EACH LINK * 

* GIVEN THE MATERIAL AND AREA DATA FROM COMMON BLOCK? * 

* * 
**«*************«**»»»*«*»»»*****•*****»•••*•**»*•****»***•*••»•*•••»»»• 

REALMS A(2),D(2),THETA(2),ALPHA(2),AREA,E,G. 

+ IYY,IZZ,MU(2),GRAV,S(4,4),P 

COMMON /BLOCKl/ A, D, THETA, ALPHA, AREA, E.G, I YY, IZZ, MU, GRAV 

DO 10 1=1,4 

10 



END 



SUBROUTINE VECTR (X,R) 

* * 

* THIS SUBROUTINE CALCULATES THE POSITION VECTOR R * 

* * 
************************************************************************ 

REAL*8 A ( 2 ) , D( 2 ) , THETA ( 2 ) , ALPHA ( 2 ) , AREA, E , G , 

+ IYY,IZZ,MU(2),GRAV, 

+ R(4,1),L,X 

COMMON /BLOCKl/ A, D, THETA, ALPHA, AREA, E,G, I YY, IZZ, MU,GRAV 
L»A(2) 

* CREATING THE LOCAL R POSITION VECTOR 

DO 10 J=l,4 
10 R(J,1)=0,D0 

R(l,l) « l.DO 
R(2,l) = X 

END 



SUBROUTINE COEFPl (X,P1) 

A*********************************************************************** 
A A 

* THIS SUBROUTINE WILL CALCULATE THE VALUES FOR THE PHIl COEFFICIENT * 

A A 

aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 



REALMS A ( 2 ) , D ( 2 ) , THETA ( 2 ) , ALPHA ( 2 ) , AREA , E, G , 

+ IYY,IZZ,MU(2),GRAV, 

+ P1(4,6),L,X,K 

COMMON /BLOCKl/ A, D, THETA, ALPHA, AREA, E,G, I YY, IZZ, MU, GRAV 

L=A(2) 

K = X/L 
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DO 10 J=l,4 
DO 10 M=l,6 

P1(J,M)*0.D0 
Pl(2,l) = 1. + K 

Pl(3,2) - -2.*K**3 - 3.*K**2 f 1. 
Pl(4,3) « Pl(3,2) 

Pl(4,5) « L*K*(K**2 + 2.*K + 1.) 

Pl(3,6) » Pl(4,5) 

END 



SUBROUTINE BINT ( TERM, XL, XU, 2, Y, NR, NC, INDEX, DD, MAT) 



GIVEN THE TEMPORARY MATRICES 2, AND Y THE INDEX # * 
CORRESPONDING TO INTEGRAL DEFINITIONS, THE DIMENSIONS OF THE • 
RESULTANT MATRIX PRODUCT (ROW, COL), AND THE H FUNCTION INDEX (D) * 
THIS SUBROUTINE WILL INTEGRATE THE RESULTANT PRODUCT MATRIX * 
OVER THE LENGTH OF THE LINK USING GAUSSIAN QUADRATURE METHODS * 
AN ASSUMPTION OF CONSTANT WEIGHT DENSITY AND CROSS SECTIONAL • 
AREA HAS BEEN MADE WRT THE LINK/ ACTUATOR VOLUME * INTEGRAL. * 



INTEGER NR, NC, INDEX, TERM 

REALMS A(2),D( 2),THETA(2),ALPHA( 2) ,AREA,E,G, 

IYY,I22,MU( 2),GRAV, 

^ AA,B,C,XL,XU,2(NR,NC) , Y(NR,NC),MAT(NR,NC),DD(4,4) 

COMMON /BLOCKl/ A, D, THETA, ALPHA, AREA, E, G , I YY, 122, MU, GRAV 
AA=.5D0*(XU+XL) 

B=XU-XL 

C*. 480 144928248768 12D0*B 

CALL INTPRD ( TERM, INDEX, DO, AA+-C, NR, NC, Y) 

CALL INTPRD ( TERM, INDEX, DD, AA-C, NR, NC, 2 ) 

DO 10 1=1, NR 
DO 10 J=1,NC 

10 MAT(I,J)=.50614268145188130D-1*( Y(I, J)+2(I, J) ) 

C=. 3983332387068 1337D0*B 

CALL INTPRD ( TERM, INDEX, DD, AA+C, NR , NC, Y ) 

CALL INTPRD ( TERM, INDEX, DD, AA-C, NR, NC, 2 ) 

DO 20 1=1, NR 
DO 20 J=1,NC 

20 MAT(I,J)=MAT(I,J)+. 11119051722668724D0*( Y(I, J)>2(I, J) ) 

C=. 26276620495816449D0»B 
CALL INTPRD ( TERM , INDEX, DD, AA+C, NR , NC, Y ) 

CALL INTPRD ( TERM, INDEX, DD, AA-C, NR, NC, 2 ) 

DO 30 1=1, NR 
DO 30 J=I,NC 

30 MAT(I, J)=MAT(I, J)+. 15685 332293894364D0 * ( Y ( I , J ) +2 ( I , J ) ) 

C=. 9171732124782490D-1*B 
CALL INTPRD ( TERM, INDEX, DD, AA+C, NR , NC, Y ) 

CALL INTPRD ( TERM, INDEX, DD, AA-C, NR , NC, 2 ) 

DO 40 1=1, NR 
DO 40 J=1,NC 

40 MAT(I, J)=B»(MAT(I,J) + . 18134189168918099D0*( Y(I, J)+2(I, J) ) ) 

DO 50 1=1, NR 

DO 50 J=1,NC 

IFdNDEX.EQ. DTHEN 

MAT(I,J)= MAT(I,J) 

ELSE 

MAT(I,J)=(MU(1)*AREA) * MAT(I,J) 

ENDIF 

50 CONTINUE 
RETURN 
END 
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SUBROUTINE DFORM( H,K,UTERM^ DC) 

************************************************************************ 
* * 

* THIS PROGRAM CALCULATES THE DEFORMATION TRANSFORMATION MATRIX * 

* FOR THE KINETIC AND POTENTIAL ENERGY CONTRIBUTION OF POINT LOAD * 

* * 
************************************************************************ 

REALMS UTERM(6),DC(4,4) 

DO 10 I»l,4 
DO 10 J=l,4 

10 DC(I,J)«0.D0 



IF(M.EQ» DTHEN 
DC(1,1)=1.D0 
DC(2,2)=1.D0 
DC(3,3)=^1.D0 
DC(4,4)*1.D0 
ENDIF 

IF((K.EQ.l) .OR. (K.EQ.7) )THEN 
DC(2,1)=UTERM(1) 

ENDIF 

IF((K.EQ.2) .OR. (K.EQ.7) )THEN 
DC(3,1)=UTERM(2) 

ENDIF 

IF((K.EQ.3) .OR. (K.EQ.7) )THEN 
DC(4,1)«UTERM(3) 

ENDIF 

IF((K.EQ.4) .OR. (K.EQ.7) )THEN 
DC(3,4)=-UTERM(4) 
DC(4,3)=UTERM(4) 

ENDIF 

IF((K.EQ.5) .OR. (K.EQ.7) )THEN 
DC(2,4)*«UTERM(5) 

DC(4,2)= UTERM(5) 

• ENDIF 

IF((K.EQ.6) .OR. (K.EQ.7) )THEN 
DC(2,3)*-UTERM(6) 
DC(3,2)=*UTERM(6) 

ENDIF 

RETURN 

END 



SUBROUTINE CCOEF(EPSLON, IOTA,H) 

************************************************************************ 

* FORMULATING THE COEFFICIENTS USED IN THE SEQUENTIAL INTEGRATION * 

* ADDITIONAL DEFINITION OF VARIABLES 

* Al-7 THE COEFFICIENTS FOR NUMERICAL INTEGRATION 

************************************************************************ 

REALMS EPSLON,IOTA,H, 

4- A0,A1,A2,A3,A4,A5,A6,A7 

COMMON /BLOCKS/ AO , A1 , A2, A3 , A4 , A5 , A6 , A7 
AO » l.DO/EPSLON/(H*H) 

A1 =IOTA/EPSLON/H 
A2 = l.DO/EPSLON/H 
A3 = 0.5/EPSLON - l.DO 
A4 = lOTA/EPSLON - l.DO 
A5 = H/2.D0 * (A4 - l.DO) 

A6 = H*(1.D0 - IOTA) 

A7 * IOTA * H 
END 
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SUBROUTINE TRANSP ( TERM, Ml, R, C, MIT ) 

* 

* GIVEN A MATRIX H THIS SUBROUTINE MILL RETURN THE TRANSPOSITION 

• OF THE MATRIX, MT 



*****^***********-k****m’tt**mtimitit’ttititinntimnitin,in, 



A*********************** 



INTEGER R,C,TERM 
REAL-8 M1(R,C), M1T(C,R) 



* 

* 

* 

A 

A 



DO 10 J=1,C 
DO 10 K=*1,R 

M1T(J,K)=M1(K, J) 

10 CONTINUE 
END 



SUBROUTINE RINT ( I . H. NOME, J, K) 

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
A A 

* THIS SUBROUTINE WILL PRINT ANT OF THE ARRAYS USED, IT IS * 

* PRIMARILY USED FOR CHECKING THE RETURNED VALUES. * 

A A 

aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 



INTEGER I,J,X,L,N,0 
REALMS M(I,J,K) 

CHARACTER'S NAME, NOME 

NAME=» • 

NAME=NOME 
DO 10 L » 1,1 

WRITE(6,1000) 'ARRAY ',NAME,' FOR I = *,L 
DO 10 N « 1,J 

WRITE (6, 1100) (M(L,N,0),0*1,X) 

10 CONTINUE 

1000 FORMAT ( IX, A8 , : AS , A9 , I 1 ) 

1100 FORMAT (IX, 7(: 2X,E9.3)) 

* END 



SUBROUTINE REFORM 

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 



* REFORMULATING THE COEFFICIENT MATRICES OF THE SMALL MOTION * 

* EQUATION OF MOTION FOR SOLUTION BY SEQUENTIAL INTEGRATION * 

A A 

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 



REAL*8 MQQ,MQN(S) ,MNN(S,S),GN(S,S),XN(S,S),FN(6),FQ, 

+ MN ( S , S ) , FNN ( 6 ) , FQTEMP , MQNTMP ( 6 ) , TEMP 1(6,6),PRD(S,S) 

COMMON /BLOCKS/ MQQ, MQN, MNN, GN, XN, FN, FQ 
COMMON /BLOCK?/ MN,FNN 

FQTEMP=( 1 . DO /MQQ ) *FQ 
DO 10 1 = 1, S 

10 FNNC I ) =FN( I ) -MQN( I ) • FQTEMP 

DO 20 1=1, S 

MQNTMP ( I ) = ( 1 . DO / MQQ ) *MQN ( I ) 

20 CONTINUE 

CALL MATMUL(PRD,MQN,6, 1, MQNTMP, 1,S, TEMPI) 

DO 30 1=1, S 

DO 30 J=1,S 

30 MN(I,J)=MNN(I,J)~TEMP1(I,J) 



END 
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SUBROUTINE TRANSF(TIME) 



THIS SUBROUTINE CREATES THE HARTENBERG-DENAVIT MATRIX FOR EACH 
JOINT AND THE TRANSFORMATION MATRICES RELATING EACH OF THE JOINT 
BODY FIXED COORDINATE SYSTEMS 

ADDITIONAL DEFINITION OF VARIABLES 
HD THE HARTENBERG DENAVIT TRANSFORMATION MATRIX 






k********************* 



REALMS A ( 2 ) , D ( 2 ) ^ THETA ( 2 ) , ALPHA ( 2 ) , AREA , E, G, 

^ IYY,IZZ^MU(2) ,GRAV,HD( 2,4,4) ,Q,QDOT,Q2DOT,TIME 

COMMON /BLOCKl/ A,D,THETA, ALPHA, AREA, E,G, lYY, IZZ, MU, GRAV 
COMMON /BLK15/ Q,QDOT,Q2DOT 
COMMON /BLOCK2/ HD 

* CREATING THE HARTENBERG-DENAVIT MATRIX FOR EACH I 
DO 10 1=1,4 
DO 10 J=l,4 



HD( 1,I,J)=0.D0 
10 HD(2,I,J)=0.D0 



HD(1,1,1) 


a 


l.DO 


HD(1,2,2) 


B 


COS(Q) 


HD(1,2,4) 




SIN(Q) 


HD(1,3,2) 


a 


SIN(Q) 


HD( 1,3,4) 


a 


-COS(Q) ■ 


HD(1,4,3) 


a 


l.DO 


HD(2,1,1) 


a 


1. 


HD(2,2,1) 


a 


A(2)»COS(THETA(2) ) 


HD(2,2,2) 


= 


COS(THETA(2) ) 


HD(2,2,3) 


= 


-SIN(THETA(2) ) 


HD(2,3, 1) 


a 


A(2)*SIN(THETA(2) ) 


HD(2,3,2) 


= 


SIN(THETA(2) ) 


HD(2,3,3) 


= 


COS(THETA(2) ) 


HD(2,4,4) 


■ 


l.DO 



RETURN 

END 



SUBROUTINE MATMUL ( PRD, MATl, R1 ,C1, MAT2, R2, C2 , PROD ) 



THIS PROGRAM RECEIVES TWO MATRICES TO BE MULTIPLIED TOGETHER 
ENSURES PROPER DIMENSIONALITY AND RETURNS THE PRODUCT MATRIX 




INTEGER R1,C1,R2,C2 

REALMS MAT1(R1,C1),MAT2(R2,C2),PR0D(R1,C2),PRD(R1,C2),SUM 

• CHECKING THE DIMENSIONALITY OF THE MATRICES TO BE MULTIPLIED 

IF (C1.NE.R2) PRINT*, 'ERROR IN MULTIPLICTION DIMENSIONS' 

* MULTIPLYING THE MATRICES TOGETHER 

DO 30 J=1,R1 
DO 20 K=1,C2 
SUM=0.0 
DO 10 H=1,C1 



IF(ABS(MAT1(J,M) ) .LE. l.E-30) MATl ( J , M ) =0 . DO 





CONTINUE 

IF(ABS(MAT2(M,K) ) .LE. 


l.E-30) MAT2(M,K)=0.D0 




CONTINUE 

SUM = SUM + MAT1(J,M) 


* HAT2(M,K) 


10 


CONTINUE 






PRD(J,K)=SUM 




20 


CONTINUE 




30 


CONTINUE 




50 


DO 50 1=1, R1 
DO 50 J=1,C2 

PROD(I,J)=PRD(I,J) 





END 
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SUBROUTINE WMATRX 



THIS SUBROUTINE MULTIPLIES THE HARTENBERG-DENAVIT MATRICES 
TOGETHER TO OBTAIN THE COORDINATE TRANSFORMATION MATRIX 

ADDITIONAL DEFINITION OF VARIABLES 



HD THE HARTENBERG DENAVIT TRANSFORMATION MATRIX 

V TRANSFORMATION MATRIX 



REALMS W(4,4),HD(2,4,4), tfSUM ( 4 , 4 ) , SUM 



COMMON /BLOCX2/ HD 
COMMON /BLOCX3/ W 



10 



DO 10 L=l,4 

DO 10 M=»l,4 

WSUM(L,M)=0.D0 
tf (L,M)=0.D0 



DO 40 J=l,4 
DO 30 X=»l,4 
SUM=0.0 
DO 20 M»l,4 

SUM « SUM + HD(1,J,M) * HD(2,M,X) 
20 CONTINUE 

WSUM(J,X)»SUM 
30 CONTINUE 

40 CONTINUE 



DO 50 I=»L,4 
DO 50 J=l,4 

50 tf(I, J)=WSUM(I,J) 

END 
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APPENDIX E 



FIGURES 



LFMR SCAN GEOMETRY 




Figure 1 N-ROSS Terrestrial Track 
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NROSS/LFMR 

BASELINE CONFIGURATION (DEPLOYED) 




Figure 2 LFMR Antenna Subsystem Configuration 
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SPACECRAFT 

BOOM 




ROTATIONAL POSITION (DGRS) 



BOOM POSITION CONSTANT TORQUE 



ZERO TIP MASS 




Figure 3 The Large Motion Position 0, Constant Torque, 
Zero Tip Mass. 
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BOOM POSITION CONSTANT TORQUE 



44 KILOGRAM TIP MASS 




Figure 4 The Large Motion Position 0, Constant Torque, 
44 Kilogram Tip Mass. 



76 



BOOM VELOCITY CONSTANT TORQUE 

ZERO TIP MASS " 



Figure 5 




(s/ava) AxiooiaA ivNoixvioa 



The Large Motion Velocity 0 Dot, Constant 
Torque, Zero Tip Mass. 
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TIME (S) 




(s/avH) AiDoiaA ivMoiivioa 



Figure 6 The Motion Velocity 9 Dot, Constant Torque, 
44 Kilogram Tip Mass. 
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TIME (S) 



BOOM ACC. CONSTANT TORQUE 

ZERO TIP MASS 




,OT*(s/avH) NoiivaaiaoDV avNOUvxoa 



Figure 7 The Large Motion Acceleration, 9 Double Dot, 
Constant Torque, Zero Tip Mass. 
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TIME (S) 



o 

Ci 




0*098 0*0 0*098- 



(2**S/aVH) 3DV IVMOIIVXOH 



Figure 8 



The Large Motion Acceleration 0 Double Dot, 
Constant Torque, 44 Kilogram Tip Hass. 
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TIME (S) 




TIP DEFORMATION CONSTANT TORQUE 

44 Kir.OGRAM TIP MASS 




,_OT. (w) NoisNaixa ivixv 



Figure 9 The Boom Tip Deflection in the Axial Direction, 
Constant Torque, 44 Kilogram Tip Mass. 
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TIME (S) 



0 

01 




,0T* (w) iiMawaovidsia a avooa 



Figure 10 The Boom Tip Deflection in the Y Transverse 
Direction, Constant Torque, 44 Kilogram Tip 
Mass . 
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TIP DEFORMATION CONSTANT TORQUE 

44 KILOGRAM TIP MASS 




O'S O'O O'S- 



j.oT* (K) iNaiMaovidsia z ivdot 



Figure 11 The Boom Tip Deflection in the Z Transverse 

Direction, Constant Torque, 44 Kilogram Point 
Mass . 
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TIME (S) 





Figure 12 The Boom Tip Torsional Deflection, Constant 
Torque, 44 Kilogram Tip Hass. 
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TIME (S) 




TIP DEFORMATION CONSTANT TORQUE 

44 KILOGRAM TIP MASS 




,.0T* (Qva) Moiivioa sixv a ivooi 



Figure 13 The Boom Tip Rotation About the Local Z 

Direction, Constant Torque, 44 Kilogram Tip 
Mass . 
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TIME (S) 



o 




Figure 14 The Boom Tip Rotation About the Local Y 

Direction, Constant Torque, 44 Kilogram Point 
Mass . 
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TIME (S) 



LEVEL I 



LEVEL_1I 



LEVEL III 





TRANSF 

HD 




WDERIV 

V WDOT W2DOT 
£(W’’W) 



COEFPl 



co EFqi . 



VECTR 



MATMUL 



TRAHSP 



6 IH .T. 

integrate 
Boom Vol. 



INTPRD 

provide BINT 
kernel 



Figure 15 Simulation Code Flow Chart 
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